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Preface 


The purpose of this book is to provide an introduction to and an 
overview of numerical modelling of the ocean and the atmosphere. It 
has evolved from a course given at Stockholm University since 1997 
and is intended to serve as a textbook for students in meteorology 
and oceanography at the master level. A prerequisite is a background 
in mathematics, physics, some geophysical fluid dynamics and pro- 
gramming. Focus will be on numerical schemes for the most common 
equations in oceanography and meteorology as well as on the stability, 
precision and other basic numerical properties of these schemes. We 
will use as simple equations as possible that still capture the properties 
of the primitive equations used in the general circulation models. For 
simplicity, the equations will often be referred to as the hydrodynamic 
equations since the numerical methods to be described here are valid 
for modelling both the ocean and the atmosphere. Due to the non- 
linearity of these equations, it is not possible to find analytical solu- 
tions. The equations are therefore instead solved numerically on a grid 
by discretisation, and the derivatives of the differential equations are 
replaced by finite-difference approximations. This is what constitutes 
a numerical model, which often is referred to as a general circulation 
model when it represents the three-dimensional (3D) global circulation 
of the atmosphere or the ocean. These models are based on the Navier- 
Stokes equations (including the Coriolis effect) and a tracer equation 
for heat in both the atmosphere and ocean and tracer equations for 
humidity and salt in the atmosphere and ocean, respectively. A cou- 
pled atmospheric and oceanic general circulation model represents the 
core part of an Earth System climate model. 

The focus here will be on the basic numerical methods used for 
oceanographic and atmospheric modelling. For more detailed and com- 
prehensive books see e.g. Durran (2010) for the numerical methods 
and e.g. Kalnay (2003) for data assimilation and numerical weather 
prediction. 
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The book starts with a short summary of the historical background 
of the numerical modelling of the ocean and atmosphere. In Chapter 2, 
we present the most common types of Partial Differential Equations 
(PDEs) that occur in meteorology and oceanography and how they can 
be classified. Chapter 3 introduces some of the most commonly used 
finite differences in both time and space and how accurate they are 
compared to the derivative. In Chapter 4, we show how to undertake a 
stability analysis using, for simplicity, the advection equation. The sta- 
bility is shown to depend both on which finite-difference schemes that 
are used and on the advection speed, time step and grid size. Some 
schemes turn out to be unconditionally unstable, but can still be used 
as a first time step. In Chapter 5 we show how, in addition to the phys- 
ical mode, a numerical mode arises and how this mode can be damped 
or suppressed. The accuracy of the numerical phase speed is examined 
in Chapter 6 and how this depends on the wave number and grid reso- 
lution. The shallow-water equations are discretised in Chapter 7 using 
three different spatial grids. Many limitations, which previously have 
been studied for the advection equation, are here examined again, but 
for the shallow-water equations. In Chapter 8 we investigate the dis- 
cretisation of the friction and the diffusion terms and how this affects 
the stability of the solution. Iterative methods for the Poisson and 
Laplace equations are demonstrated in Chapter 9. Implicit and semi- 
implicit schemes are shown in Chapter 10. Chapter 11 deals with the 
semi-Lagrangian technique for the advection equation. Different model 
coordinates for atmospheric as well as oceanic models are presented in 
Chapter 12. Using a highly simplified approach, 3D-modelling is intro- 
duced in Chapter 13. Chapter 14 gives a brief description of how some 
atmospheric general circulation models use spectral methods as “hor- 
izontal coordinates”. Some “pen-and-paper” theoretical exercises are 
given in Chapter 15 and Chapter 16 is devoted to a number of GFD 
computer exercises. 

We wish to thank Laurent Brodeau, Bror Jönsson, Joakim Kjellsson, 
Gurvan Madec as well as our students for valuable input during all these 
years. 


Kristofer Döös, professor of climate modelling, 

Peter Lundberg, professor emeritus of physical oceanography 
and Dr. Aitor Aldama Campino, all at the Department of 
Meteorology, Stockholm University 


1. Introduction 


1.1 What is a numerical model of the circulation of the 
atmosphere or the ocean? 


A numerical model of the circulation of the atmosphere and/or the 
ocean is basically constituted by 


A grid covering the spatial domain under consideration. 
e Discretised equations describing the conservation of 
momentum, mass, energy and salt /moisture. 
e Open and/or solid boundaries. 
e A specified initial state of the system. 


The schematic cartoon shown in Figure 1.1 illustrates a system of this 
type pertaining to a coupled ocean-atmosphere model with correspond- 
ing forcing and the exchange between the two components of the model. 


1.2 Brief historical background 


A short summary is here given of the historical background to the 
numerical modelling of the ocean and atmosphere. 

Bjerknes (1904) was the first to discuss the possibility of predicting 
and modelling the circulation of the atmosphere. For this purpose he 
proposed seven equations with seven unknown variables: 


e Three equations of conservation of momentum for the three 
velocity components based on Newton’s second law. 
e The continuity equation, i.e. the conservation of mass. 
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Figure 1.1. Schematic illustration of a numerical ocean-atmosphere 
circulation model 
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e The equation of state for ideal gases. 

e The equation of conservation of energy based on the first law of 
thermodynamics. 

e A conservation equation for water mass in the atmosphere. 


These equations have become known as the “primitive equations” since 
they do not deal with a filtered quasi-geostrophic equation but go back 
to Bjerknes’ original formulation of the problem. These equations also 
require boundary conditions at the top and bottom of the atmosphere 
as well as the initial state of the atmosphere. The problem is that 
although we have an equal number of equations and variables, as well 
as appropriate boundary and initial conditions, these equations cannot 
be solved analytically. 

The first attempt to solve the equations numerically was made by 
Lewis Fry Richardson in his spare time when serving as an ambulance 
driver during the First World War. Richardson (1922) applied finite- 
differencing techniques to the equations suggested by Bjerknes. 

A decade earlier Richardson (1911) had been the first to propose 
the use of finite differences when approximating partial differential 
equations, insights he put to use when attempting numerical weather 
prediction. His first numerical weather prediction (NWP) failed, how- 
ever, due to two reasons. First, no computers were available and 
Richardson estimated that it would require a staff of 64,000 persons 
to compute a 24-hour forecast in 24 hours. The second reason, which 
was revealed by Lynch (2006), was due to a failure to apply smoothing 
techniques to the data, which led to his calculations yielding erroneous 
results. 

These numerical instabilities were independently investigated 
by Courant, Friedrichs and Lewy (Courant et al., 1928), pure 
mathematicians who examined techniques of solving general partial 
differential equations using finite differencing. They found certain con- 
ditions regarding the choice of time step and grid size that must 
be respected for the numerical solution to be stable. These results 
were later further developed by Charney, Fjørtoft and von Neumann 
(Charney et al., 1950), the trio that in the late 1940s made the first 
successful numerical weather forecast, based on integration of the con- 
servation equation for the absolute vorticity. 
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This single-equation approach of following the motion of an air col- 
umn instead of using the seven-equation set proposed by Bjerknes 
(1904) was suggested by Carl-Gustaf Rossby, who was also the driv- 
ing force behind the first real-time numerical weather prediction made 
in two runs 23-24 March 1954 in Stockholm by Harold Bedient and 
Bo Döös (Döös and Eaton, 1957; Wiin-Nielsen, 1991; Persson, 2005a). 

The rationale underlying the need of numerical weather predic- 
tion is immediately obvious; less so when ocean forecasting is con- 
cerned, where two severe inundations of the northwestern European 
coast served as the driving agent. In the Netherlands, severe flooding 
in 1916 with around 300 fatalities led the government to entrust the 
physicist and Nobel Laureate Hendrik Lorentz with the task of devel- 
oping predictive techniques. This ultimately led to Lorentz formulat- 
ing a method for storm-surge forecasting in the 1920s. An even worse 
instance of North-Sea coastal flooding took place in 1953, with around 
2000 deaths, which led to the Hamburg oceanographer Hansen (1956) 
developing a numerical storm-surge model based on finite differenc- 
ing of the shallow-water equations including the non-linear advection 
terms. Hansen focused on the fast barotropic gravity waves instead of 
filtering them out as was done in numerical weather-prediction mod- 
els. Subsequently another German oceanographer (Fischer, 1959) con- 
structed a numerical shallow-water model for the North Sea with the 
finite-difference schemes described in detail. 

The return to the full primitive equations for NWP models, which 
Richardson (1922) had used in the very beginning, was inevitable since 
the quasi-geostrophic equations, although very useful for understand- 
ing the large-scale extratropical dynamics of the atmosphere, were not 
accurate enough to allow continued progress in NWP. A complete 
Atmospheric General Circulation Model (AGCM) based on the primi- 
tive equations was developed by Smagorinsky (1963) at the Geophys- 
ical Fluid Dynamics Laboratory (GFDL) in Princeton. Here also a 
large part of the subsequent model development took place, including 
the first Ocean General Circulation Model (OGCM), developed only a 
few years later by Bryan and Cox (1967). Syukuro Manabe and Kirk 
Bryan combined their models to yield the first coupled Atmosphere- 
Ocean General Circulation Model (AOGCM), cf. Manabe and Bryan 
(1969), which they used for climate studies. 

Ocean and atmosphere GCMs have since then increased in number 
all over the world and have improved in many aspects; higher reso- 
lution using more powerful computers, better parameterisations, more 
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observations to feed the models with. The basic numerics have also 
improved, but the fundamental properties of the finite differences and 
the numerical methods remain. This is why focus in this book will be 
on the limitations of the numerical schemes that to a large extent were 
discovered during the 20th century. 

For a more comprehensive view of the historical background, there 
are number of studies focusing on the evolution of NWP, e.g. Persson 
(2005a,b,c) and Wiin-Nielsen (1991). 


2. Partial Differential Equations 


Before going into the details of numerical procedures, we first classify 
the most usual types of Partial Differential Equations (PDEs) that 
occur in meteorology and oceanography. Partial differential equations 
are of vast importance in applied mathematics and engineering since so 
many real physical processess can be modelled by them. Second-order 
linear PDEs of the type needed for modelling the ocean and atmosphere 
circulation can be classified into three categories: elliptic, parabolic and 
hyperbolic. These partial differential equations in two dimensions are 
general linear equations of second order: 
ou ou ou ðu o 


u 
H H H = 0. 2.1 
ana Pavia oye Oe E fu+g=0 (2.1) 


This equation has two independent variables (x and y) and one depen- 
dent (u) as well as two unspecified functions (f and g). We will show 
here that by transforming x, y and u into new variables, it is possible to 
apply Equation 2.1 to any second-order linear PDE and to classify it. 

As the reader may note, this expression bears somewhat of a 
resemblance to the equation for a conic section: 


ax? + bxy + cy? + dz + ey + f =0, (2:2) 


where a, b, c, d, e and f are constants. As illustrated by Figure 2.1, 
this algebraic equation represents an ellipse, parabola or a hyperbola 
depending on whether b? — 4ac is negative, equal to zero or positive, 
respectively. 

This indicates that one analogously can classify the PDE according 
to Table 2.1. 
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Figure 2.1. A diagram showing conic sections, viz. an ellipse, a parabola 
and a hyperbola. Note that these geometrical representations of Equation 
(2.2) are not solutions to the general PDE given by Equation (2.1). 


Table 2.1. Classification of the three PDE types and their different kinds of 
boundary conditions. A Dirichlet boundary condition specifies the function 
on the boundary of the domain. A Neumann boundary condition specifies 
the value of the normal derivative of the function on the boundary of the 
domain. A Cauchy boundary condition specifies the values of the function 
and its normal derivative on the boundary of the domain. 


PDE b? — 4ac Boundary & initial conditions Examples 
Elliptic | b? — 4ac < 0 | Dirichlet /Neumann/Cauchy Laplace Eq. 


Parabolic | b? — 4ac = 0 | One initial + two boundary Diffusion Eq. 
conditions 


Hyperbolic | b? — 4ac > 0 | Two initial + two boundary Wave Eq. 
conditions 


Each type of system is associated with significantly different char- 
acteristic behaviour, and the solution scheme for each type of equation 
can also differ. All three classes of PDEs are represented among the 
most common equations in hydrodynamics and require the specification 


Partial Differential Equations 9 


of different kinds of boundary conditions. We will now provide some 
examples of these equations. 


2.1 Elliptic equations 


The two-dimensional (2D) versions of the Laplace and Poisson 
equations are classical examples of elliptic equations, representing e.g. 
the relationship between the stream function and vorticity or the 
steady-state temperature in a plate: 


Pu Ou 


Jz T ay 0 (or g (u, x,y). (2.3) 


In conformity with Equation (2.1) it is recognised that a = 1, b = 0, 
c= ] and d = e = f = 0, which leads to b? — 4ac = —4 < 0, and hence 
the PDE must be elliptic. 

Another example of an elliptic equation in hydrodynamics can be 
found from the simplest equations of motion in the atmosphere and 
the ocean, viz. the linearised shallow-water equations without friction. 
These equations can be both hyperbolic (as shown in Section 2.3) and 
elliptic when in steady state as will be shown here. 


ðu an 

Ov On 

—, z= — — 2.4 
ore By’ (2.4b) 
Oh ðu Ov 


Here u and v are the velocity components in the x and y direction, 
respectively, 7 is the surface-height perturbation from the undisturbed 
depth H of the fluid, g is the acceleration of gravity and t represents 
the time. Above f = 2Q sin ọ is the Coriolis acceleration, where Q is 
the angular frequency of the Earth’s rotation and y the latitude. In 
what follows f is set to be a constant. 

Note that the variables x and y in Equation (2.1) have been desig- 
nated in an arbitrary fashion and do not have anything to do with the 
independent variables x and y in Equations (2.4). From these equations 
it can be deduced that 


ð | & o? o? 
z E | f? — gH (= + =) | n= 0. (2.5) 
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If we integrate this equation in time we obtain 


82 oO? 82 
a E f’n on (so | Z) n = § (x,y), (2.6) 


where (x,y) is an integration constant independent of t. When the 

problem is stationary, the geostrophic relationship is obtained: 

0? o? 
+ 

Ox? 0y? 


Peon ( ) n = £ (x,y). (2.7) 
In analogy with Equation (2.1), we find that a = —gH, b = 0 and 
c = —g H, which leads to b? — 4ac = —4g° H? < 0, and hence this PDE 
must be elliptic. 


There are three different types of possible boundary conditions for 
this class of PDEs: 


e u specified on the boundary (Dirichlet), 

e 0Ou/O7, where 7 is the normal vector, specified on the boundary 
(Neumann), 

e au+ 0u/O7 specified on the boundary (Cauchy). 


2.2 Parabolic equations 


An example of a parabolic PDE is the diffusion equation for the depen- 
dent variable T corresponding to e.g. temperature, salinity, humidity 
or any passive tracer: 


oT 0 
Ot ðr? 
Based on Equation (2.1), a = K, b = 0 and c = 0, leading to b? — 4ac = 


0, and hence the PDE is parabolic. 
Assume that T (x,t) is the temperature distribution along the x-axis 


K >00. (2.8) 


of a conducting rod as a function of time. To solve the equation over 
the interval 0 < « < L one must prescribe an initial condition T (x, 0) 
on 0 < a < L. The boundary conditions T (0, t) and T (L, t) must also 
be specified during the whole time period under consideration. 


2.3 Hyperbolic equations 


This class of PDEs describes wave motion. A typical example is the 
equation governing the vibrating string or the gravity waves in the 
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ocean or atmosphere. In Section 2.1, we found that the shallow-water 
equations were elliptic in their steady state form. Here we will show 
that they can be hyperbolic in their time-dependent form. The wave 
equation can be derived from the simplest possible set-up of the 
shallow-water equations: 


Ou On 

a eg 2: 
On Ou 

dl ray pea 2. 

Ot Ox ea 


By eliminating u between the two equations we obtain the wave 
equation: 


n g 


ae 7" aa" 
In analogy with Equation (2.1) it is found that b? — 4ac = 4gH > 0 
and hence this equation is a hyperbolic PDE. 


(2.10) 


The following first-order PDE, known as the advection equation, 
can also be classified as hyperbolic, since its solutions satisfy the wave 
equation: 

a S (2.11) 
where co = yvgH is the speed of gravity waves. By first differentiating 
the advection equation with respect to t and then a second time with 
respect to x, it is possible to eliminate 8’u/ðtðx between the resulting 
two equations and we obtain the wave equation. 

In order to obtain unique solutions of Equation (2.10) for 0 < xz < 
L we need boundary conditions at x = 0 and x = L, viz. the ends 
of the spatial domain, and two initial conditions. Possible boundary 
conditions are u(0,t), u(L,t), Ou(0,t) /Ox or Ou(L,t) /Ox specified 
at x = 0 and x = L. Necessary initial conditions are u(z,0) and 
Ou (x, 0) /Ot specified over 0 < x < L. 


2.4 Overview 


Having introduced these three types of PDEs, it is important to under- 
line that the behaviour of the solutions, the proper initial and/or 
boundary conditions, and the numerical methods that can be used 
to find the solutions depend essentially on the type of PDE dealt 
with. Although non-linear multidimensional PDEs in general cannot 
be reduced to these canonical forms, we need to study these prototypes 
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in order to develop an understanding of their properties, and then 
apply similar methods to the more complicated equations governing 
the motion of the ocean and atmosphere. 


Exercises: 


1. Solve the following equation by separation of variables: 


a (2.12) 
dt ax’ l 
with the initial condition 
u (x,0) = —Ae”**. (2.13) 
2. Show that the advection equation 

Ou Ou 
— = -—c— 2.14 
dt Ox ee) 


has the general d’Alembert solution u = f(x — ct), where f is 
an arbitrary continuously differentiable function. Interpret the 
equation geometrically in the xt-plane. Solve the equation with 
the following initial condition: u (x,0) = g(x). 

3. Derive the wave equation from the shallow-water equation 


system 
av 
a 2.1 
Bt gVN, (2.15a) 
On > 
— =-HV. V 2.1 
ot i (215b) 


where g is the gravitational acceleration, 7 the free-surface 
deviation from the equilibrium depth H, V the horizontal 
velocity vector (u,v) and V = (0/0z,0/Oy). 


3. Finite Differences 


3.1 The grid-point method 


Let us study a function u with one independent variable x: 
w= u(x). 


Suppose we have an interval L, which is partitioned by N + 1 equally 
spaced grid points (including the two at the limits of the interval). 
The grid length is then Az = L/N and the grid points are located at 
x; = jAx, where j = 0,1,2,...N are integers. Let the value of u at x; 
be represented by uj. 


3.2 Finite-difference schemes 
The formal mathematical definition of the derivative of a function 
u(x) is 
du, = u(a@+ Az) — u(x) 
eect Ar 


, (3.1) 


which is illustrated in Figure 3.1. A derivative and its finite-difference 
approximation differ in that in the latter case Ax remains finite and will 
not tend to zero. Note that in what follows, the term “finite difference” 
will be taken as synonymous with finite-difference approximations of 
derivatives. 

We will now derive expressions which can be used to give an approx- 
imate value of a derivative at a grid point in terms of grid-point values. 
The finite-difference schemes can be constructed between values of u; 
over the grid length Az. As illustrated in Figure 3.2, the first derivative 
of u (x) can be approximated in three ways: 
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x x+Ax 
Figure 3.1. As Ax decreases, the red secant approaches the blue tangent, 


which is the derivative of the function u (green) at the point x when 
Aa —> 0. A finite difference is when Ax Æ 0. 


. Äx , (j-1)Ax  jAx (j+1)Ax 


X 


Figure 3.2. The Euler-backward, -forward and centred finite-difference 
approximations of the derivative of a function u (x) at the point jAx. The 
function is defined at the grid points x = jAx so that uj = u (x) = u (jAz), 
where Az is the grid length and j = 0, 1,2, ... are integers. Note that we 
have chosen an exceptionally difficult point around which to approximate 
the derivative, this to highlight potential difficulties with the method. 


forward-difference scheme: 


du E! Uj+1 — Uj 
dz), ~ Ar ’ 
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centred-difference scheme: 


du\  Uj+1 — Uj- 
dx } , ~ 2Az ° 


backward-difference scheme: 


du a Uj — Uj-1 
dz), ~ Ar ` 


These various schemes introduce errors that can estimated by deriv- 
ing the finite-difference schemes in a more rigorous way using Taylor 
expansions. The Taylor series for f (y) around y = a is 


f(y) =F (a)+ly — a) F (a)+5 (y — a)? f” (a)+ + (y — a)” F (a). 


2 
(3.2) 
Note that a Taylor series is an expansion around a specific point, which 
yields a function that can be evaluated at another point. 


3.2.1 Forward-difference scheme 


Using the Taylor expansion given by Equation (3.2) and substituting 
f(y) by u(x), a by z; and y by x41, it is found that 


du 1 2 (du 1 3 (du 
Wee ery es pais z Gi = ay 
Uj4+1 = Uz + (te) +ga) (cs) te (= 3 


The forward difference can now be expressed as 


Uj41 — Uj du 1 du 1 2 (du 
= Ax) | — —(A — 
Ar (aa z) (= fe! a) dz? } , 


1 dtu 1 du 
Shay le) bea ee, E 
tor?) (Z) +a Ae) (a) of) 


4 


The difference between this expression and the real derivative 
(du/dx), is 


1 du 1 2 [Ëu 
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which is denoted the truncation error associated with the approxima- 


et 


tion of the derivative. The terms that have been truncated, i.e. “cut 
off”, are represented by £. Hence we have an accuracy of first order 


with 

e = O (Az), (3.6) 
which is the lowest order of accuracy that is acceptable. Note that 
we have assumed that Az is small so that (Ax)™ << (Ax)” where 


m is greater than n. The higher-order terms in Ax can therefore be 
neglected. 


3.2.2 Centred-difference scheme 


The accuracy of the centred-difference scheme can be determined from 
Equation (3.3) and the analogous representation of uj—1: 


7 du 1 o {du 1 3 [u 
Uj-1 = Uj Ax (T) +500" (Ta) 5am (= , 


so that 


Uj41 — Uz—1 du 1 2 ( du 1 4 {du 
= A — (A — reia 
2Ax (ue) ta a G ee) \ aes), 


Here the truncation error is of second order: 


5 : (Az)? (=) $2. = 0 [(An)'], (3.8) 


3.2.3 Centred fourth-order difference scheme 


A scheme with fourth-order accuracy can be obtained if we undertake 
the Taylor expansion given by Equation (3.2) and substitute f(y) by 
u(x), a by x; and y by 254: 


du 1 of du 1 3 ( du 
es 2 z a = (2A nals 
Upto = Uj + 2AT (Ge) +5 2a) ean x) (= 


1 4 [d'u 1 s [ču 
= CII de AGP fe aw. 
Ta CD) (=) Pa C42) | ae ; We 


J 
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and analogously, 


du 1 2 [Pu 1 a [Ëu 
,-2 = Uj — 2Ax | — = (2A — ] —-(2A — 
es E. a ey ol 2) (= j 


Uj+2 — Uj—2 du 4 2 Bu 16 4 du 
= poker | a a 
4Ar (a) g (A) (= TIl D ‘a 


This scheme is, as the previous centred scheme, accurate to second 
order, and if we combine the two centred schemes so that 


A Uj41 — Uj- 1 Uj42 — Uj- = (=) 1 (Ax)* (=) COPEN 
J j 


3 2Azx 3 4Az 
we find an accuracy of fourth order, i.e. € = O [(Az)’]. 


3.2.4 Centred-difference scheme on a staggered grid 


Finite differences on staggered grids (cf. Chapter 7) often require that 
the approximations be located between two contiguous grid points: 


=) Uj41 — Uj 
- py nT, (3.12) 
(z Jti Ax 


The accuracy of this scheme is obtained by first making a Taylor expan- 
sion around #;41/2 and then taking into account the point 7;4,, which 
leads to 


n Ar (=) mn 1 (5) (5) 
Uj+1 = Uj+1/2 aN See 
' 2 \dz/ i 2\2 dz? J yaya 


1 (Az \?® (a 
| ( z) (5 3 Eaa D 
6\2 dx pais 


and then using the same Taylor expansion but with regard to the point 


x; so that 
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Ax (=) 
Uj = Uj41/2 — F 2 
7 3 2 dz) jiy 


1 (Ax? (du 1 /Axr\? (Bu 
+ ( ) ( ) ( ) ( ) ta. (3.14) 
2\ 2 da?) a 6\2 da?) j 


Subtracting the latter equation from the former and dividing by Ax 


leads to 
se = u d 1 (Az\? (& 
ce a U ( +) 4 ( z) ( =) +.. (3.15) 
Az dx faye 6\ 2 dx janie 
Hence the truncation error is of second order: 


e= : (4°) ar +... = O [(Az)’]. (3.16) 


Note that the only difference between this truncation error and the one 
for the previous centred finite difference given by Equation (3.8) is that 
this one is reduced, which is simply due to that the distance between 
the points is Ax instead of 2Az. 


3.3 Time-difference schemes 


The time-derivative schemes that are used for PDEs are often rel- 
atively simple, usually of second-order and sometimes even only of 
first-order accuracy. There are several reasons for this. First, it is a 
general experience that schemes constructed to have a high order of 
accuracy in time are mostly not very useful when solving PDEs. This 
is in contrast to the experience with ordinary differential equations, 
where very accurate methods, such as the Runge-Kutta scheme, are 
extremely successful. There is a basic reason for this discrepancy. With 
an ordinary differential equation of first order, the equation and a single 
initial condition is all that is required for an exact solution. Thus, the 
numerical-solution error is entirely due the degree of inadequacy of the 
scheme. With a PDE, the error associated with the numerical solution 
arises from both the shortcomings of the scheme and the insufficient 
information about the initial conditions, which only are known at dis- 
crete grid points. Thus, an increase in accuracy of the applied scheme 
improves only one of these two components, and the result is not too 
impressive. 
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Another reason for not using a finite-difference scheme of high 
accuracy of the time derivatives is that, in order to meet a stability 
requirement of the type to be discussed in the next chapter, it is usu- 
ally necessary to choose a time step significantly smaller than that 
required for adequate accuracy. Once a time step has been specified, 
other errors, e.g. from the spatial differencing, are much greater that 
those due to the time differencing. Thus, computational effort is better 
spent in reducing these errors, and not in increasing the accuracy of 
the time-differencing schemes. This, of course, does not mean that it 
is not necessary to carefully consider the properties of various possible 
time-differencing schemes. Accuracy is only one important considera- 
tion when choosing a scheme. 

To define some schemes, we consider a general first-order differential 
equation: 
where typically u = u (x,t). The independent variables x and t are 
space and time. f is thus a function of u, x and t, corresponding e.g. 
to the advection equation where f = —cOu/0x. 

In order to discretise the equation we divide the time axis into seg- 
ments of equal length At. The approximated value of u(t) at time 
t = nAt is denoted u”. In order to compute u"t! we need to know at 
least u” and often also u”~'. A number of time-difference schemes are 
available. 


3.3.1 Two-level schemes 


These are schemes that use two different time levels, n and n + 1, so 
that the time integration yields 


(n+1)At 
ut =u” +f f (u, t) dt. (3.18) 
nAt 


The problem now is that f only exists as discrete values f” and f+! 
at times nAt and (n + 1)At, respectively. 


Euler-forward scheme 
This is defined as 


u"! =u" + At f”. (3.19) 
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Euler-forward time integration 


Figure 3.3. Time integration with the Euler-forward scheme. 


Here the truncation error is O (At), i.e. the scheme is accurate to first 
order. It is said to be uncentred, since the “time derivative” pertains to 
time level n + 1/2 and the function to time level n. As in the previous 
section this Euler-forward scheme is of first-order accuracy (Figure 3.3). 


Euler-backward scheme 
This is defined as 


yor = u” + At ft. (3.20) 


The Euler-backward scheme is uncentred in time and is accurate to 
O (At). If, as here, a value of f is taken at time level n+1 and f depends 


n+1 


on u, ie. u”**, the scheme is said to be implicit. For an ordinary 


differential equation, it may be a straightforward matter to solve for 
unt, 

For a PDE it will, however, require solving a set of simultane- 
ous equations, with one equation for each of the grid points of the 
computation region. If no value of f depends on u”+! on the right- 
hand side of the equation above, the scheme is said to be explicit. 
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In the very simple cases when f only depends on t, such as 
e.g. du/dt = —yu, the discretised equation becomes u”*t = u” + 
At (—yu"*"), which can be rearranged as u”! = u” /(1 + yAt) so that 
there are no terms at time level n + 1 on the right-hand side. In this 
case the discretised equation can be integrated despite being implicit. 


Crank-Nicolson scheme 


The Crank-Nicolson scheme is based on the trapezoidal rule. If we 
approximate f by its average between time levels n and n+1 we obtain 


1 
ut! <u" + SAt (f+ fr). (3.21) 


This scheme is implicit, since it requires information from the future 
time level n + 1. Note here that the finite-difference approximation is 
centred at n + 1/2 between the two time steps n and n + 1, viz. 


du) tP urtu” 
— = ——___— . 3.22 
(5) Ai Pe 


The truncation error £ can be found from the Taylor series in Equation 
(3.2). Substituting f(y) by u(t), a by tt"? and y by t”t!, we obtain 
the following expression: 


utt = yr te 2 du -_ atl 1 [AtN\? (Pay 
7 2 dt 2 2 dt2 
1 (At 3 Bu n+1/2 
+3($) (=) +... (3.23) 


Substituting f(y) by u(t), a by t”+!/2 and y by t”, we obtain the Taylor 
expansion of the function u (t) at time level n + 1/2: 


n n+1/2 At du cee 
u =U So 
2 \ dt 


1 (At\? (@u\"? 1 (AEN? [Bu 
| a fee (B24) 
2\ 2 dt? 6 \ 2 dt? 
Subtracting Equation (3.24) from Equation (3.23) and dividing by At 
we obtain 


n+l _ zm n+1/2 it d n+1/2 
siat -($) a (an? ( =) Fas B25) 


At dt 24 dt? 


22 Basic Numerical Methods in Meteorology and Oceanography 


which can also be expressed as : 


du n+1/2 yrtl — u” 
— — 2 
@ ae te (3.26) 
where 
1 Guy 
pem (At)? (>) +...= O [(At)’]. (3.27) 


The truncation error £ is therefore of second order for the Crank- 
Nicolson scheme. 


Matsuno’s forward-backward scheme 


To increase the accuracy compared to the Euler-forward and -backward 
schemes we can construct iterative schemes such as the Matsuno 
scheme, which is initiated by an Euler-forward time step: 


unt =u" + At f”. (3.28) 


In this case the value of the obtained u”*! serves as an approximation 
of f"*1, which hereafter is used to make a backward step to yield a 
final u”! : 


got = u” + At f, (3.29) 
where 


fE = f (unt, (n + 1)At). (3.30) 


* ? 


This scheme is explicit and of accuracy O (At). 


Heun scheme 


This is similar to the Matsuno scheme and is explicit but of second- 
order accuracy, since the second step is made using the Crank-Nicolson 
scheme: 


uT! =u" + Atf”, (3.31) 


At 
ut ut Sf + ith). (3.32) 


The mid-point scheme 

The mid-point scheme, also known as a second-order Runge-Kutta 
method, is like the Matsuno and Heun schemes a multi-stage scheme, 
which uses an intermediate estimate of the solution throughout the 
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time-step. The scheme consists of first taking an Euler-forward scheme 
one half time step: 


A 
urtt/? = yr + — f” (3.33) 


and then using the solution at the mid-point time level n + 1/2 to 
integrate with a centred scheme: 


u” = u” + At fte, (3.34) 


The fourth-order Runge-Kutta scheme 

The fourth-order Runge-Kutta scheme is similar to the previous scheme 
but integrates the solution with four steps instead of two. The scheme 
consists of first applying an Euler-forward scheme one half time step: 


At 
uu E fh, (3.35) 


and hereafter recomputing the solution at time level n + 1/2 with an 


Euler-backward scheme making use of u?*+1/?: 


At 
ut =u” + T oz (3.36) 


Subsequently a centred scheme is used to arrive at time level n + 1 
making use of u"*1/?: 


ge = u” + At f+, (3.37) 


The final solution is obtained by making use of Simpson’s rule with the 
four values u”, u?tt/2, u”+t1/2 and u?t!, resulting in 


At 
ut! = u” + a a Ope ag ne 4 ee (3.38) 


This fourth-order Runge-Kutta method is stable, and as its name 
indicates, accurate to fourth order. It is seldom used except in e.g. 
some regional high-resolution NWP models. Runge-Kutta methods of 
even higher order exist, but are rarely used in ocean or atmosphere 
circulation models. 


3.3.2 Three-level schemes 


These schemes use the time at three levels and the time integration 
becomes 


(n+1)At 
yet = u”! + f f(u,t)dt. (3.39) 
( 


n—1)At 
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The simplest three-level scheme is to assign f a constant value equal 
to that at the middle of the time interval of length 2At, which yields 
the leap-frog scheme 


unt! = ut" + 2At f” (3.40) 


of accuracy order (At)*. This has been a widely used scheme in both 
atmospheric and oceanic models. Many of today’s atmospheric cir- 
culation models use, however, Lagrangian time-stepping, which will 
be introduced in Chapter 11. Fourth-order schemes are, as mentioned 
above, sometimes also used. In some models such as the ROMS ocean 
model, there are several different time-schemes that can be used. 


Exercise: 


Determine the order of accuracy of the centred discretisations of the 
advection equation 


4. Numerical Stability 


Here we will study partial differential equations (PDEs) with one 
dependent and two independent variables. Intuitively one often thinks 
of the advection equation as describing the evolution of a tracer drifting 
passively with the flow. However, we shall here consider various sim- 
plified forms of the advection equation describing the advection of a 
dependent variable. In practice, this has proved to be the most impor- 
tant part of the hydrodynamic equations for the atmosphere and the 
ocean. 

We will use the advection equation to investigate what is required 
of the numerical schemes to yield stable solutions to the PDE, i.e. 
solutions such that small perturbations do not grow in time, but rather 
decrease. 


4.1 The advection equation 


This hyperbolic advection equation (cf. Section 2.3) is 


Ou Ou o 
Ot Ox f 
where u = u(x,t) and c is the prescribed phase speed. Analytical 


solutions are of the form u(x,t) = ue", where c = w/k, w being 
the frequency and k the wave number. 

Let us now consider one among many possible discretisations of this 
equation by using a centred difference scheme in both time and space: 


n+1 n-1 n aj? 
j ly pestle, Wel 


Uj j | 
2Az 


T = U. 4.1 
2At j cee) 
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0 1 j-i j jn ea 
Figure 4.1. A finite-difference grid in time and space. The grid points used 
by the centred schemes in both time and space are shown as green dots. The 
red dots are the necessary initial-condition points and the blue dots 


illustrate the two boundaries, at which conditions must be prescribed. The 
brown dot shows the “here and now” point 7,7. 


Note that all the grid points for this scheme, which are shown by the 
green dots in Figure 4.1, are located around the “here and now” point 
j, n shown as the brown dot in the graph. When used in model appli- 
cations, this needs to be reformulated in order to integrate values to 
the “future” time step n + 1: 


oo =a =p (uta _ u?) f (4.2) 
where 
cAt 
= 4.3 
u m (4.3) 


is known as the Courant number, sometimes denoted the CFL-number 
(Courant et al., 1928). Figure 4.2 shows how a cosine “hump” propa- 
gates, when integrated 25 time steps using Equation (4.2). 
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4.2 Initial and boundary conditions 


We also need initial and boundary conditions in order to integrate 
Equation (4.2) forward in time on a grid such as that shown in 
Figure 4.1. The initial condition is that all the u values (red dots) 
have to be prescribed, i.e. the variable u must have specified values 
at time step n = 0. Furthermore we need to be able to integrate a 
first time step, which is not possible with a three-level time scheme. 
We therefore use an Euler-forward scheme for the first time step and 
hereafter proceed using the leap-frog scheme for the rest of the time 
integration. 

The boundary conditions imply that values for all the uj and u7, 
(blue dots in Figure 4.1) have to be prescribed, i.e. the variable u must 
have values at all time steps on the two boundaries located at j = 0 
and j = JX, where JX is the total number of grid cells. 

When we have a periodic domain, where e.g. the values on the 
eastern boundary equal those on the western boundary, which is the 
case for global models of the Earth, we can use periodic boundary 
conditions. For practical purposes this boundary can e.g. be located 
at the Greenwich meridian where the longitude can be expressed 
as 0° or 360°. When u is computed in e.g. a Fortran code, at 
the eastern (j=0) and western (j=JX) boundaries one will need 
u-values for “j=-1” and “j=JX+1”. This can easily be accomplished 
in the 7-loops of the model code by introducing jm=j-1 and when 
jm=-1, replacing it by JX-1. The same procedure is used for the 
eastern boundary with jp=j+1 and when jp=JX+1, it is replaced 
by jp=1. 

A segment of a Fortran code of the numerical time integration of the 
discretised advection equation (4.2) could look like this: 


mu=1. ! Courant number 

u(:,0)=1. ! initialise the field to something, e.g. 1 
! time loop 

do n=1,NT-1 


do j=0,JX 
jp=j+t 
if (jp==IJX+1) jp=1 
jm=j-1 
if(jm==-1) jm=JX-1 


u(j,n+1) = u(j,n-1) - mu * (u(jp,n)-u(jm,n)) 
enddo 
enddo 
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Figure 4.2. The advection equation integrated numerically with centred 
schemes in both time and space using Equation (4.1). The first time step 
has been integrated with an Euler-forward scheme. The initial value is 
u?=° = cos [m (j — 50)/20] for 40 < j < 60 and u?™ = 0 for the remaining 
grid points. The Courant number for this scheme is u = cAt/Az = 1. 
Solutions are graphed for time steps n=0, 5, 10, 15, 20, 25. Note the 
“numerical noise” propagating in the “wrong” direction in the form of 
small-amplitude jagged waves, a topic to be further considered in 
Chapter 6. 


4.3 Stability analysis of the leap-frog scheme 


In order to study the stability we will use what is known as the von 
Neumann method, which in fact was first briefly introduced by Crank 
and Nicolson (1947) and later more rigorously by Charney et al. (1950), 
in which latter study von Neumann actually was the last author. This 
method is generally not possible to use for non-linear equations and one 
is therefore limited to apply it to linearised versions of the equations 
in a numerical model. In general a solution of a linear equation can be 
expressed as a Fourier series, where each Fourier component is a solu- 
tion. We can thus test the stability using solely one Fourier component 
of the form 


u? _ ge Ae-Epn At) 
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where k is the wave number. Note that the phase speed Cp, with the 
index D pertaining to finite differencing, is an approximation of the 
phase speed c occurring in the differential equation. It is this Cp, result- 
ing from solving Equation (4.1), which here will be investigated. We 
note that 


ik|jAe—Cp(n+1)At] — ,n,-ikOpAt _ „n 
=urye =u, 


n+l 
Uj 


= Uge 


—ikCpât 


where A= e is the amplification factor. Similarly, for n — 1 we 


find that 
wre (4.4) 


These results can be generalised as 
UnT™ = up A™ (4.5) 


and 


ur = Auges, (4.6) 


From this we can deduce that if |u| is not going to “blow up” when 
integrating in time, it is required that 


[Alle hee | < 1 (4.7) 


and reversely, if |A| > 1 the solution is unstable and “blows up”. For 
the stability condition to be fulfilled, Cp must be real. This technique 
for determining stability, based on examining the amplification factor 
A, is known as the von Neumann method. We also need expressions for 
the Fourier components of the spatial derivatives: 


ik|(j+1)Ae—CpnAt] _ grae 4.8 


Sg ety (4.9) 


J 


nm 
Uji = Uge 
n — ik|(g-1)Aw—CpnAt 
Uji = Uge [ ] 


Let us now return to Equation (4.1) and introduce A: 


n — \—1,,n ikâAzy n _ p—ikār,;mn 
Aur AÀ u; | P uy —e ur o, (4.10) 
2At 2AT 


which, since e** — e~** = 2isina, can be simplified to the quadratic 


equation 


X? + 2iusin (kAr)\—-1=0. (4.11) 
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Furthermore x?+az + 8 = 0 yields that x = —a/2+ ,/a?/4 — 6, and 
hence Equation (4.11) has the solution 


\ = —ipsin (kA) + y- (usin (kAz)}? + 1. (4.12) 


Keeping in mind that the absolute value of a complex number a + ib is 


Va? + b?, we find that if 
[usin (kâr)? <1, (4.13) 


then 


|A|? = [usin (kKAx)]? + \v- [usin (kAx)]? + i =, (4.14) 


i.e. the scheme we are presently examining is stable if Equation (4.13) 
holds true, a condition which also can be formulated as 


| usin (kAvr)| < 1. (4.15) 
Since |sin (kAzx)| < 1, we have conditional stability if 


Ax 


<1 < 
|u| S Lor |e s S, 


(4.16) 


this since we require stability for all wave numbers k and have hence 
chosen the “worst” case, i.e. when sin (kA) = 1 . The reverse case is 
when 


[usin (kâr)? > 1. (4.17) 


For simplicity we define 
a = psin(kAz), (4.18) 
so that Equation (4.17) is reduced to 
a >l. (4.19) 


It is then preferable to rewrite Equation (4.12) as 


A= ta bie? — 17% (-a +a — 1) l (4.20) 


or equivalently 


2 
|A? = (a + Va? =T) = 20? + 2a =T- 1, (4.21) 
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which has at least one root larger than one. Consequently when pz > 1 
(c > Axv/At), then |A| > 1 for at least some wave numbers k and thus 
the solution is unstable and hence of limited practical use. 

The leap-frog scheme is thus conditionally stable. Given a spatial 
resolution Ax we require a time step At not exceeding Aa/c for 
Equation (4.16) to be valid for the fastest possible phase speed in the 
system as illustrated by Figure 4.3. This is known as the Courant- 
Friedrichs-Lewy (CFL) criterion (Courant et al., 1928). In this par- 
ticular case of the advection equation with centred finite differences 
both in time and space, the Courant number should not exceed 1 
(u = cAt/Ax < 1) in order for the CFL condition to be satisfied. 
This criterion will, as we shall see later, change depending on which 
equation and finite difference scheme we are dealing with. 


Figure 4.3. The Courant-Friedrichs-Lewy (CFL) stability criterion for 
schemes centred in both time and space. 
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An example of a successful integration of the advection equation 
with centred finite differences in space as well as in time has previously 
been given in Figure 4.2, where we see how the initially prescribed 
trigonometric “hump” with an amplitude equal to 1 progresses right- 
wards along the x-axis. This representation of a time-dependent process 
is, however, somewhat unwieldy and in practice what is known as a 
Hovmöller (1949) diagram is most frequently used, cf. Figure 4.4. 
This highly compact visualisation is based on graphing the the spatio- 
temporal evolution of the process in a coordinate system with the 
ordinate representing time and the abscissa pertaining to the spatial 
evolution of the process. The left-hand panel of Figure 4.4 thus shows 
precisely the same results as those given in Figure 4.2, where the grad- 
ually weakening coloration of the fringes of the “diagonal bands” do 
perfect justice to the spatial structure of the trigonometric “hump” in 
Figure 4.2. The right-hand panel of Figure 4.4 shows how for a Courant 
number larger than 1 instability sets in after less than 10 time steps 
and subsequently grows to encompass the entire spatial range under 
consideration. 
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Figure 4.4. Hovmöller diagrams of the advection equation integrated 
numerically with a leap-frog scheme in time and a centred scheme in space, 
viz. Equation (4.1). The first time step has been integrated with an 
Euler-forward scheme. The initial value is u?=° = cos [m(j — 50)/20] for 

40 < j < 60 and uno = 0 for the remaining grid points. The Courant 
number is u = cAt/Az = 1 in the left panel, which yields a stable solution, 
but in the right panel the solution is seen to “blow up” for u = 1.1. Note 


that periodic boundary conditions have been applied. 
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Here it may finally be noted that if the advection equation dealt with 
here using centred finite differencing is applied to the ocean, then we, 
for a given spatial resolution Ax, must adjust the time step At to satisfy 
u < 1. Here the relevant phase speed is the one for long gravity waves 
c = /gH, where H is the depth of the ocean. In order to guarantee 
numerical stability we need to determine the maximum depth Hmax 
and set the the time step to satisfy At < Ar/c = Ar//gHyyax- 


4.4 Euler-forward scheme in time 


Let us now apply an ad hoc numerical scheme to the advection 
equation, so that instead of Equation (4.1) we have 


n+1 n n n 
Uz -— u; Uj+ı — Uj-1 
H = 0, 4.22 
At TA 2A : (Pa 
which yields 
à= 1 — isin (kAz). (4.23) 


As before the stability criterion is that |A| < 1, and since the absolute 
value of X satisfies 


|A|? = 1+ [usin (kAx)]’ , (4.24) 


it is recognised that |A| > 1. The solution consequently grows with 
time, independently of how we choose the time step. This scheme is 
thus unconditionally unstable, which is simply referred to as unstable. 


4.5 The upstream scheme 


This scheme is denoted upstream or upwind since it looks for infor- 
mation from where the wind or current comes by using two different 
spatial schemes depending on the direction of the velocity. 


u — u” u? — u? 
1 T L e~ Ap =0 ife>d, (4.25a) 
n+l n n n 
Uj — uj Uji — Uj ; 
= 0. 4.25b 
Ai eS ae 0 ifc< ( ) 


If c > 0 one should use the backward scheme in space from Equation 
(4.25a) in order to use information from where the flow comes. 
Reversely if c < 0 then one should use the forward scheme in 
space given by Equation (4.25b). Figures 4.6 and 4.7 show the time 
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Figure 4.5. The advection equation integrated with an Euler-forward 
scheme in time and a centred scheme in space. In the left panel the Courant 
number is u = 0.9 and in the right panel 0.4 . Irrespective of the choice of 
Courant number, the solution will sooner or later “blow up”. 


evolution when Equation (4.25) is integrated with different Courant 
numbers u. 
A von Neumann stability analysis yields 


A =1-p[l1-—cos(kAz) + isin (kAz)], 
resulting in 


A? = {1 — u [1 — cos (kAx)]¥ + [usin (kAzr)? 
= 1 +2u [|1 — cos (kAzx)| (u — 1). (4.26) 


It is unfortunately not straightforward to see directly from this equation 
the corresponding stability criterion. For this we need to graph its 
amplification factor À as a function of both resolution kAzx and Courant 
number u as in Figure 4.8. This shows that |A|? <lfor0O<yp<1l, 
which also implies that c must be positive to ensure stability. The 
scheme is hence conditionally stable. This is similar to the leap-frog 
scheme but with the major difference that the amplification factor A 
will be reduced when u < 1, which leads to a decrease of the amplitude 
at every new time step. This deamplification is clearly visible in the 
upper right panel of Figure 4.6 as well as in the right panel of Figure 4.7, 
where the upstream scheme has been integrated with = 0.5. The 
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Figure 4.6. The advection equation integrated numerically with the 


upstream scheme given by Equation (4.25a). Stable solutions are shown in 
the upper panels, where to the left u = 1 and to the right u = 0.5, the latter 
solution being clearly dissipative. The lower panels show unstable solutions 
with to the left u = 1.1, and to the right u = —1, where in the latter case 
the spatial forward-difference scheme given by Equation (4.25b) should 
instead have been used. 


decrease in amplitude is stronger for the short waves as illustrated 
by Figure 4.8. The upwind scheme should hence be integrated with 
a large Courant number u, i.e. as close to one as possible. In this 
very simple case of the advection equation it is possible to use u = 1, 
but in a more comprehensive model, such as a GCM, this will not be 
possible and hence the scheme will lead to dissipation. The upwind 
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Figure 4.7. The amplitude u} for different time steps n, computed with the 
upstream scheme given by Equation (4.25a), with u = 1 to the left and with 
u = 0.5 to the right. Note the rapid deamplification when integrated with 

u = 0.5. 


Figure 4.8. The squared amplification factor of Equation (4.26) as a 
function of the Courant number yp for different wavelengths measured in the 
number of grid lengths Ax. Note that there is a minimum at u = 0.5, which 
leads to a deamplification of, in particular, the short waves. 


advection scheme was proposed by Courant et al. (1952) and used in 
early numerical weather-prediction models due to its good stability 
properties. It still finds use in idealised oceanic box models as well as 
in some general circulation models. The finite-volume model developed 
at ECMWF is e.g. using the upwind scheme in its advection-transport 
algorithm making use of an iterative approach (Smolarkiewicz et al., 
2016). When using an upwind scheme one should, however, realise 
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that it is highly diffusive. It uses backward spatial differencing if the 
velocity is in the positive x-direction, and forward spatial differencing 
for negative velocities. The term upwind denotes the use of upwind, 
or upstream, information when determining the form of the finite 
difference scheme; downstream information is ignored. 


4.6 Stability analysis of the fourth-order scheme 


Let us now study an advection scheme that uses the spatial scheme 
accurate to fourth order given by Equation (3.11): 


n+1 n—1 n n n n 
uj =u; 4 Uj — Uj-  1Uj+— Uj-2 
= 0. 4.27 

za S (; 2Ax 3 4Az aoa 


A von Neumann stability analysis yields 
Nee if [8 sin (kAx) — sin (2kAz)] À — 1 = 0, (4.28a) 
which has the solution 


A=- if [8 sin (kAx) — sin (2kAzx)} 


2 


ž y 2 [£ [8 sin (kAz) — sin (2kAz)]| +1. (4.28b) 


If 
l [8 sin (kAx) — sin (2kAx)]] <i, (4.29) 


then 


2 


(a= [$ [8 sin (kAx) — sin (2kA2)]| 


+ fy- [£ [8 sin (kAzx) — sin (2kAz)]| Í + y =1, (4.30) 


i.e. this scheme is stable if Equation (4.29) is fulfilled. This condition 
can also be expressed as 


£ [8 sin (kAz) — sin (2kAx)] <1 


or 
r 6 
F S 8sin (kAx) — sin (2kAx) 


(4.31) 
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Figure 4.9. The right-hand side of Equation (4.31) as a function of kAa/n. 


The minimum value in red, which is the stability criterion with a Courant 
number p = cAt/Arx ~ 0.729. 
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Figure 4.10. The advection equation integrated with the fourth-order 
spatial numerical scheme given by Equation (4.27). The solution is stable in 
the left panel where u = 0.728, but unstable in the right panel when 

u = 0.730. 
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The exact stability criterion is found by determining the kAzx for which 
the right-hand side of Equation (4.31) attains its minimum as illus- 
trated by Figure 4.9. This turns out to be for kAz/r = 0.57, which 
leads to the scheme being stable for the Courant number u < 0.729. 
Figure 4.10 shows how the integration with the fourth-order scheme 
evolves in time with a Courant number just above and below its criti- 
cal value. 


Exercises: 


1. Consider the leap-frog scheme for the advection equation: 


n+l n-1 n n 
a Uj41 — Uj- 
=0. 4.32 
Ar COKE (4.82) 
Use 
u? = A per (4.33) 


and show that the amplification factor is 


A = —ipsin (kAx) + yı — ( usin (kAz)}’. (4.34) 


2. Show that for u > 1 in the previous exercise, we will have one 
of the solutions to the differential equation “blowing up”, at 
least for some wavelengths. 

3. Discretise the advection equation with Euler-forward schemes 
in both time and space. Show that for c > 0 (backward 
scheme), the amplitude of the solutions will grow in time (these 
thus being unstable). But for c < 0 (forward scheme) the 
amplitude will decrease in time, i.e. the solution is stable. 

4. Undertake a stability analysis of the following discretisation of 
the advection equation: 


a) (4.35) 


Is it centred or uncentred? 


5. The Computational Mode 


Here we shall examine some of the consequences of partial differential 
equations having been approximated with finite-difference schemes. 


5.1 The three-level scheme 


One problem with a three-level scheme such as the leap-frog one is that 
it requires more than one initial condition to start the numerical inte- 
gration. From a purely physical standpoint a single initial condition for 
u”™= should suffice. However, in addition to this physical initial condi- 
tion, for computational purposes three-level schemes require an initial 
condition also for w"=' (in what follows the index j will be omitted 
since there is no spatial dependence in the equations to be dealt with). 
This value cannot be calculated using a three-level scheme, and will 
generally have to be determined using some type of two-level scheme. 


Consider the oscillation equation: 


= = wu, (5.1) 
where u = u(t). The analytical solution is 
u = ue. (5.2) 
Equation (5.1) can be integrated numerically with a leap-frog scheme: 
utt =u! + 2iwAtu". (5.3) 
If we now examine the amplification factor (cf. Chapter 4) we find 


A? — 2iwAtA — 1 =0, 
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which has the solution 


Ai2 = iwAt+ \/ — (wAt)? +1. 


Thus there are two solutions of the form u”t! = Au”. Since we are 
dealing with a linear equation, its solution will be a linear combination 
of these two: 


up =Afuy ; uy = Agua, 

so that 

u” = adPul + bAP us, (5.4) 
where a and b are constants. If u”! = Au” should represent the approx- 
imation of the true solution, then A > 1 when At —> 0. This condition 
is satisfied by A,, but also for A. — —1. The solution involving ), 
is denoted the physical mode and that with A> the computational or 
numerical mode induced by the numerical scheme. This latter mode 
changes sign for each even and odd n. 


A straightforward way to illustrate the computational mode is to 
study the simple case when w = 0, viz. 


du 

— =0, 

dt 
which has the exact solution 

u(t) = constant. 
Discretisation using the leap-frog scheme yields 
ues aia (5.5) 

For a given physical initial condition u"=° = C, we consider two 


special choices of u"=?: 


1. If calculating u”™t happens to yield the true value C4, then for 
all n 


utt = u” (5.6) 
or 
unt? = dou”. 


In this case we have obtained a numerical solution identical to 
the true solution and which consists solely of the physical mode. 
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2. Suppose now instead that when calculating w"=! we obtain 
u=! = —u”™., Then for all n 
ut = —y" (5.7) 
or 


amt! = dow”. 


The numerical solution now consists entirely of the 
computational mode. 


5.1.1 The computational initial condition 


A suitable choice of the computational initial condition is of vital 
importance for obtaining a satisfactory numerical solution for short 
simulations, where the initial condition is crucial (as is the case for 
weather forecasts but less so for long climate simulations). The com- 
putational initial condition (u"“'), which is one time step ahead of the 
physical condition (u”~°), can be calculated with a single Euler-forward 
time step. Although this scheme is computationally unstable, it can be 
used for a single time step since a considerable number of steps are 
required before the solution finally “blows up”. 

The computational initial condition for our academic oscillation- 
equation case is then, with an Euler-forward time step: 


u =u? + iwAtu”™. (5.8) 


An alternative computational initial condition, useful when the solu- 
tion is less sensitive to the initial condition, is to assign the same value 
to both time steps (u"=! = u"=°), but, as shown above, this can imme- 
diately trigger a computational mode. 


5.2 Suppression of the computational mode 


The computational mode induced by the leap-frog scheme can be sup- 
pressed in two different ways: 


5.2.1 Euler-forward or -backward schemes at regular intervals 


The computational-mode problem can be solved by integrating with an 
Euler-forward or -backward scheme at regular intervals constituted by 


a certain number of time steps (50 or so are often used). For Equation 


(5.5) from the previous example this would imply that u"t! = u” is 
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Figure 5.1. Illustration of the physical mode in blue from Equation (5.6) 
and the computational mode in red from Equation (5.7). The initial 
physical condition is set to u”=° = 1 and the computational initial condition 
is set to u"=! = —1. The resulting solution in black will thus change sign 
every time step since the two modes are uncoupled. 


taken every 50 time steps, which eliminates the computational mode 
so that in Equation (5.4) a= 1 and b = 0. 


5.2.2 The Robert-Asselin filter 


Another way of suppressing the computational mode, which is the most 
common one used in atmospheric models, is to employ a Robert-Asselin 
filter (Robert, 1966; Asselin, 1972). First one applies a leap-frog inte- 
gration to 


oY = g(u), (5.9) 


this in order to obtain the solution at time level n + 1: 
urtt = ue? + 2Atg(u"), (5.10) 


whereafter the filter is applied as a smoothing between the three time 
levels n — 1, n and n+ 1 so that 


ue =u +y (ug) — 2u” +u), (5.11) 


where the index f indicates the filtered values and y is the Asselin 
coefficient, usually chosen to range between 0.01 and 0.2. The next 
“frog jump” will be 


ut? =a Arg (ut). (5.12) 
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A segment of a Fortran code of the numerical time integration of the 
discretised advection equation (4.2) employing a Robert-Asselin filter 
could look like this: 


gamma = 0.01 ! Robert-Asselin filter coefficient 
mu = 1. - gamma ! Courant number 
! time loop 
do n=1,NT-1 
! Leap-frog scheme of the advection equation 
do i=1,JX-1 
u(€i,nti) = u(€i,n-1) - mu * (uCiti,n)-u(i-1,n)) 
enddo 
! Robert-Asselin filter 
do i=1,JX-1 
u€i,n) = u(i,n)+gamma*(u(i,n-1)-2.*u(i,n)+u(i,n+1) ) 
enddo 
enddo 


Note that the added term resembles smoothing in time; an approx- 
imation of an ideally time-centred smoother is 


up = u” +y (a — 2u” a (5.13) 


In our particular case of the discretised oscillation given by Equation 
(5.3) we can estimate the damping effect of the Robert-Asselin filter by 


wnAt into the smoother 


introducing the discretised solution u” = ue 
(Equation (5.13)), with the exception that u”~' is taken as an unfiltered 


value. This results in 
u% = u” [1 — 4ysin? (wAt/2)]. (5.14) 


The computational mode, the period of which is 2At, is hence reduced 
by (1 — 4y) every time step. Since the field at n — 1 is replaced by 
an already filtered field, the Robert-Asselin filter introduces a slight 
difference compared to this simplified filter. 

Another drawback of the Robert-Asselin filter is that it affects the 
stability of the schemes. A stability analysis can be performed in the 
case of the advection equation discretised with centred schemes (same 
as Equation 4.2): 
=u} — H (U — Uj) (5.15) 
with the filter of Equation (5.11). The amplification factor then 
becomes 


n+l 
Uj 


A=-itat+y+J/(1-7)? -@’, (5.16) 
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where a = psin(kAzx). For stability we require as usual that |A| < 0, 
which is obtained when u < 1—7¥ or u+y < 1. The adding of a Robert- 
Asselin filter results hence in a stricter stability condition and requires 
a shorter time step for a given spatial resolution than without filter. 
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Figure 5.2. Hovmöller diagrams of the advection equation, integrated with 
the leap-frog scheme and a Robert-Asselin filter. The solution at the first 
time step has been set to equal the initial condition instead of integrating 
with an Euler-forward in order to generate a computational mode. Top left 
panel with no filter and an oscillating solution due to the computational 
mode. Top right panel with an unstable condition (u + y < 1). Lower left 
panel with too much filtering but stable. Lower right panel with just enough 
Robert-Asselin filter to smooth the solution and a Courant number to 
match y for stability. 
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In Figure 5.2, we have illustrated the Robert-Asselin filter for 4 
different cases, where the advection equation has been integrated with 
the leap-frog scheme of Equation (5.15). The initial condition is a cosine 
wave. The first time step has not been integrated as usual with an 
Euler-forward scheme but set to be constant in time. This in order to 
immediately generate a computational mode, which we then try to filter 
out. The first case is with no filter and a persistent computational mode. 
The second case is with a filter but without respecting the stricter 
stability criterion, which makes the solution unstable and the wave 
“blows up”. The third case is with a too strong filter, which dampens 
the wave too much although it filters out the computational mode. The 
fourth case is with just enough filter to eliminate the computational 
mode without too much damping. 

Another test case is illustrated in Figure 5.3, where a shallow-water 
model has been integrated with different coefficients of the Robert- 
Asselin filter. 
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Figure 5.3. Time evolution of a shallow-water model integrated without a 
Robert-Asselin filter in black, with a Robert-Asselin filter y = 0.1 in red and 
with y = 0.01 in blue. Note that the computational mode vanishes 
completely when the stronger filter is applied, but the amplitude of the 
solution is at the same time damped. 
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5.2.3 The Robert-Asselin-Williams filter 


It was recognised that when used with the leap-frog scheme, the fea- 
ture of the Robert-Asselin filter of not conserving the mean degrades 
the numerical accuracy. Williams (2009) tackled this problem by intro- 
ducing an extra step in the filtering process in order to include the 
possibility of conserving the mean value. The resulting filter is imple- 
mented in leap-frog integrations as follows: 


ut = u% + 2Atg (ur), (5.17) 

Up = UR + (ut? — Qu? + utr), (5.18) 
b= 

upt! Sas — at= 5 a) (ut? — 2u} + uty). (5.19) 


This Robert-Asselin- Williams filter introduces an extra operation that 
is straightforward to implement and does not represent a significant 
computational expense compared to the Robert-Asselin filter. It also 
introduces a new parameter, a, such that 0 < a < 1, where a = 1 
corresponds to the traditional Robert-Asselin filter. Williams (2009) 
showed that a value of a = 0.53 minimises spurious numerical impacts 
on the physical solution and yields the closest match to the exact solu- 
tion of the equation under consideration over a broad frequency range. 


6. The Computational Phase Speed 


We shall now investigate the accuracy of the computational phase speed 
associated with using discretisations in space as well as in time. 


6.1 Dispersion due to the spatial discretisation 


Let us first examine the advection equation with a centred scheme in 
space: 


we find 
sin (kAx) 
Cp = c 
Dave ge 
where Cp is the computational phase speed and c the prescribed ana- 
lytical phase speed. Their ratio should ideally be as close as possible 
to one, but is 
Cp _ sin(kAz) 
c kde 
The computational group velocity is 


Cp, = ae ar Te (kAzx), 


which is dispersive since it depends on the wave number k, cf. 
Figure 6.1. 
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Figure 6.1. The computational phase speed associated with centred finite 
differencing in space compared to the analytical phase speed. The black line 
represnts the solution of the continuous equations which pertains the 
nondispersive analytical case, i.e. the phase speed is the same as the group 
velocity c = cg. The blue curve is the computational phase speed normalised 
by dividing with c. Note that when the wave number increases (i.e. the 
wavelength decreases) the computational phase speed deviates from the 
analytical phase speed. The phase speed is clearly dispersive since the waves 
propagate at different speeds depending on their wavelengths. The red curve 
shows the analogously normalised computational group velocity Cpg which 
at wavelengths shorter than 4 grid cells (kAx < 7/2) is in the opposite 
direction of cg. 


6.2 Dispersion due to the time discretisation 
The effects of the centred scheme in time will be analysed in the same 
way as done for the effects of the centred scheme in space: 


i —1 
Ti — u; ðu 
| 


2At 


Inserting the wave solutions 


u” (x) = ie CE ODAN, 
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we obtain 
arcsin (wAt) 
kAt i 
where Cp is the computational phase speed and c the analytical phase 


Cp = 


speed. Their ratio should ideally be as close as possible to one, but is 


Cp _ arcsin (wAt) 
c wAt ` 


The computational and analytical group velocities also differ, so that 


Co dey aaO 1 


G adk gk yA 


The computational phase speed and group velocity (both normalised 
with c) are presented in Figure 6.2, graphed as functions of wAt. 


6.3 Dispersion due to spatial and temporal resolution 


Let us now investigate the effects of the leap-frog scheme on the phase 
speed and group velocities using the advection equation with centred 


schemes in both time and space: 
n+l n-1 n nh 
Ur U a E 


2At 


By inserting a wave solution 


u? _ uget*iAe-Gpnae) 


we obtain 


1 At 
p= TE arcsin Fe sin(kAz)| ; 


where Cp is the computational phase speed and c the analytical phase 
speed. Their ratio should ideally be as close as possible to one, but is 


Cp 1 
-S arcsin |u sin (kAzx)], 6.1 
2 = -gag avesin [usin (kA) (6.1) 
where the Courant number is, as previously defined, pp = cAt/Az. 
This computational phase speed Cp is a function of the wave number 
k and the resolution Az. The finite differencing in space thus causes 
a computational dispersion. As kAzx increases, Cp decreases from c to 
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Figure 6.2. The computational phase speed associated with centred finite 
differencing in time compared to the analytical phase speed. The black line 
represents the solution of the continuous equations which is the 
nondispersive analytical case, i.e. the phase speed is the same as the group 
velocity, viz. c = cg. The blue curve is the computational phase speed 
normalised by dividing with c. Note that when the time step increases 
relative to the wave frequency w, the computational phase speed increases 
and deviates from the analytical phase speed. The red curve shows the 
computational group velocity Cpg- 


zero when kAx = 7, which corresponds to the shortest possible wave 
with a wavelength of two grid cells (A = 2Az). Thus, all waves propa- 
gate at a slower speed than the analytical phase speed c, with this decel- 
erating effect increasing as the wavelength decreases. The two-grid-cell 
wave is stationary. Note that if u = 1, which is the Courant-number 
limit of stability for the advection equation, the computational phase 
speed is the same as the analytical one, viz. Cp = c. 

The reason for the two-grid-cell wave being stationary is obvious 
when looking at the wave illustrated in Figure 6.3. For this wave we 
have uj+ı = u;_1 at all grid points, corresponding to Ou, /Ot = 0 in the 
advection equation. 
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The computational group velocity is here 


Co, = ae a ae im arcsin usin (kAc)]} : (6.2) 


A= 2Ae 


5 Sp ea AN Ilg 


Figure 6.3. The two, four, six and eight grid-interval wave with a 
wavelength of A = 2, 4,6 and 8A. 


Noting that the derivative of the inverse sine function is 


d 1 d 
pe arcsin [f (x)] = a ; (6.3) 
we obtain 
1 1 d 
Cog = [usin (kAzx)] 
At y1 — [usin (kAz)]? ah 
ccos (kAz) (6.4) 


E — [usin (kAx)]? | 


Both computational speeds are functions of the wave number, and thus 
we recognise that the spatial differencing again results in computational 
dispersion. Since the analytical group velocity is cą = c it makes sense 
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Figure 6.4. Computational dispersion of the leap-frog scheme. The curves 
show the ratio Cp/c as a function of the normalised wave number kAa/7z to 
the left and as a function of the number of grid lengths Az per wave length 
L to the right. The different curves correspond to the Courant numbers 

(u = cAt/Az) indicated on the them. The ideal solution (black line) is 
Cp/c = 1. Note that when the wave number increases (i.e. the wavelength 
decreases), the computational phase speed deviates from the analytical 
phase speed. The phase speed is clearly dispersive since the waves propagate 
at different speeds depending on their wavelengths. 


to study the ratio between the analytical and computational group 
velocities: 
Co, cos (kAx) 


Cg 7 y1 — [usin (kAx)]? 


(6.5) 


This ratio is shown in Figure 6.5. 
We have encountered two effects in this chapter. The computational 
phase speed was found to be slower than the analytical phase speed 


4 5 6 7 8 9 10 
L/IAx 


Figure 6.5. As in Figure 6.4 but for the group velocity. 
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and changes with the wave number; the wave is hence dispersive. The 
same holds true for the computational group velocity. 


Exercise: 


Derive the computational phase speed 


7. The Shallow-Water Equations 


In this chapter we will consider the equations describing the horizontal 
propagation of gravity and inertia-gravity waves. These equations are 
often referred to as the linearised shallow-water equations. We will be 
dealing with a system of two or three partial differential equations 
of first order and have two or three dependent variables (one or two 
velocities and pressure/free-surface height). The system of equations 
will always be equivalent to a single differential equation, but one of a 
higher order (which can be obtained from the system by elimination of 
dependent variables). 


7.1 The one-dimensional shallow-water equations 


We start with the simplest possible set-up of the shallow-water 
equations for a flat bottom, viz. that associated with gravity waves: 


Ou ðh 

peas eee Ll 
at Oa’ a 
ðh ðu 

— = —H— lb 
Ot Ox’ er) 


where g is the gravity, h is the thickness of the fluid and H its unper- 
turbed value, i.e. when u=0. We seek wave solutions of the form 


u (x, t) = u ™, (7.2a) 


h (x,t) = hot’, (7.2b) 
which yield the frequency equation 


w? = gHk?, 
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dh 
zZ dt 


Depth/height 


Bottom of ocean or atmosphere £ r+AT 


Figure 7.1. Schematic illustration of the 1D shallow-water equations for a 
flat bottom. 


and hence the phase speed is 


o= 7 edy 


This shows that the gravity waves can propagate along the z-axis in 
both directions with the speed gH, which is not a function of the 
wave number, and consequently the waves are non-dispersive. 


7.1.1 Spatial discretisation but continuous time derivatives 


As illustrated by Figure 7.2, there are two types of possible grids for 
these shallow-water equations. We can take the two dependent variables 
at the same points: 


Ouz ja — hj- 
a I 2r (raa) 
Oh; Uji ~ Uj-1 

= -H n 7.3b 
ot 2AT ( ) 


This is known as a unstaggered grid. It is also possible to alternate the 
grid points in space: 


Ou; Aj4a 7 hj 

= A 
t T Mr ’ vey 
Oh; Uy — Uji 

= -H ~ = 7.4b 
Ot Ax ( ) 


This is known as a staggered grid and already at this stage we can see 
that one advantage is that it reduces the number of grid points for a 
fixed truncation error. 
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Unstaggered Ax 5 


Uj-3 UY» Uji Wj Wj W42 Wj 


“j2 Hi n; Wj 
> 
x 
hi- hj Ais) 
Staggered anne enn eee > 
Ax 


Figure 7.2. Top: unstaggered grid with two dependent variables, both at 
every grid point. Bottom: staggered grid with two dependent variables at 
alternate grid points. 


The computational phase speeds and group velocities associated 
with these two schemes can be obtained by inserting wave solutions 
of the form 


uj = ripe iAe—ept) 


h; = hye’ (ih42—“pt) . 


The computational phase speed derived from the unstaggered-grid 
equations (7.3) becomes 


_ Wp, sin (kAx) 
cp = — = +ygH LAr (7.6) 


and the group velocity is 


= d(wp) d(kCo) 
Cb; = a ae +,/gH cos (kAz). (7.7) 


These are hence the same as for the spatially discretised advection 
equation in the previous chapter. The corresponding results from the 
staggered-grid equations (7.4) become 


sje Lg (7.8) 


and 


Cp; = d (wo) = OAE a) = +,/gH cos (=) ; (7.9) 
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Figure 7.3. The computational phase speed Cp associated with the centred 
finite-difference scheme in space compared to the analytical phase speed c 
for the one-dimensional shallow-water equations as a function of the 
normalised wave number kAa/z to the left and as a function of the number 
of grid lengths Az per wavelength L to the right. The black line represents 
the solution of the continuous equations which is the nondispersive 
analytical case, 7.e. the phase speed is the same as the group velocity c = cg. 
The blue curve is the computational phase speed normalised by dividing 
with c. The red curve shows the computational group velocity Cpg which at 
wavelengths shorter than 4 grid lengths (kAa < 7/2) propagates in the 
wrong direction. The green and orange curves are from the staggered grid. 


These phase speeds and group velocities are shown and inter-compared 
in Figure 7.3. Note that when the wave number increases (i.e. the 
wavelength decreases), the computational phase speed deviates from 
the analytical phase speed. The phase speed is clearly dispersive since 
the waves propagate at different speeds depending on their wavelengths. 
Apart from this the staggered grid has the advantage of reducing the 
number of grid points, and we observe that the waves with kAz > 7/2, 
which are those shorter than 4 grid lengths and have the largest phase- 
speed error, are eliminated. 


7.1.2 Spatial and temporal discretisation 


Equations (7.4) and (7.3) also need to be discretised in time in order 
to be amenable to a numerical solution. The most straightforward type 
of time differencing is the three-level leap-frog scheme. 


Unstaggered-grid case 


Discretisation of the 1D shallow-water equations using centred schemes 
in both time and space on an unstaggered grid yields 


ntl n-1 n n 
an Aa — AR 
= 7.10 
2At I- Ar Cin 
a oe. (7.10b) 


2At 7 2AT 
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We seek solutions of the form 


u? = uget tkAe—wpnat) 


h? = hoe (t Ar-upnAt) , 


which after insertion in Equations (7.10) yields the computational 


phase speed 


1 ; . 
Cp = = mma" [usin (kKAx)] , (7.12) 


where u = /gHAt/Az is the Courant number. The ratio between the 
computational and analytical phase speeds should ideally be as close 


as possible to one, but is 


Cp Cp 1 . . 
= = .1 
7 oH TUAT arcsin [u sin (kAzx)], (7.13) 


which is identical to Equation (6.1) found for the advection equation. 
The computational group velocity is also the same as for the advec- 
tion equation with the following ratio between the computational and 
analytical group velocities: 


Cbg _ cos (kAz) 
vg y1 — [usin (kAz)]? 


(7.14) 


Both computational speeds are functions of the wave number, and thus 
we recognise that the spatial differencing again results in computational 
dispersion, viz. the same result as obtained for the advection equation 
with centred schemes. 


Staggered-grid case 
The same procedure as used above, but on a staggred grid, results in 


n+l 


n-1 n n 
Us! — u; hey h} 


fj 
= Ll 
2At I Ae ’ (ela 
ae — aia Uy — Uji 
= -H J .1 
2At Ar’ (deb) 


which after inserting wave solutions yield the computational phase 
speed 


1 : , 
Cp = q = tga resin [2u sin (kAz/2)]. (7.16) 
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The ratio between the computational and analytical phase speeds 
should ideally be as close as possible to one, but is 
Cp Cp 


1 
a ee arcsin [2u sin (kA /2)] , (7.17) 


which is shown in the two lower panels of Figure 7.4 together with the 
results from the unstaggered-grid case in the upper panels. 

The ratio between the computational and analytical group velocities 
for the staggered-grid case becomes 


Cry cos (kAx/2) 


E (7.18) 
vgH y1 — [2u sin (kAa/2)]° 


which clearly differs from that obtained in the unstaggered-grid case. 
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Figure 7.4. The computational phase speed of the 1D shallow-water 
equations discretised with centred differences in time. The two upper panels 
represent the unstaggered-grid case and the lower panels the staggered-grid 
one. The curves show the ratio Cp/c as functions of the normalised wave 
number kAa/z to the left and as functions of the number of grid lengths 
Az per wavelength L to the right. The different curves correspond to the 
Courant numbers p indicated on them. 
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Figure 7.5. As in Figure 7.4 but representing the ratio of the group 
velocities from Equations (7.14) and (7.18). The two upper panels represent 
the unstaggered-grid case and the lower the staggered one. The curves show 
the ratio Cp,/c as a function of the normalised wave number kAz/z to the 
left and as a function of the number of grid lengths Ax per wavelength L to 
the right. The different curves correspond to the Courant numbers u 
indicated on them. 


Stability analysis 


The stability of the staggered-grid set of equations above can be deter- 
mined by applying a von Neumann stability analysis, which yields 


À — AT! etkan — 1 
= h .1 
2Nt Uo g Ax 03 (7 9a) 
À — )-l 1-— etkAe 


We eliminate uo and ho between these two equations and find that 


A=)" H ikAx —ikAz 
( 2At ) = ay (e^ —1) Ge ), (7.20) 


which results in two quadratic equations: 


X? + iar -—1=0, (7.21) 
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where a = 2AtV/gH/Azsin(kAx/2). The corresponding four roots 
are: 


\1,2,3,4 = tia V 1—a?. (7.22) 


The requirement for stability is that |A| < 1, which is satisfied if a < 1, 
corresponding to 


F. (7.23) 


To satisfy the stability criterion for all wavelengths, we require the 


At/gH . (=) 1 
——*— sin (| — ] < 
Az 2 


Courant number u to fulfill 


VgHAt _ 1 
= — <. 24 
u ne 5 (7.24) 


Exercise: 


Show that the stability criterion for the unstaggered-grid case is that 
the Courant number must satisfy u < 1. 


7.2 Two-dimensional shallow-water equations 


Let us now consider one of the simplest possible subsets of the equations 
of motion in the atmosphere or the ocean, viz. the linearised shallow- 
water equations (frequently denoted the inertia-gravity wave equations) 
in two dimensions: 


Ou Oh 

T fu IFz (7.25a) 
Ov Oh 

si = -g 2 
Oh ðu v 

ii E 4 z) . (7.25¢) 


Here f = 2Q sin y is the Coriolis acceleration, where Q is the angular 
frequency of the Earth’s rotation and y the latitude. In what follows f 
is set to be a constant. As before, we seek wave-type solutions: 


(u, U, h) = (uo, Vo, ho) ee, (7.26) 

Insertion into Equations (7.25) yields the following result for the 
frequency: 

w = P + gH (k? +P), (7.27) 


which describes the dispersion relationship for Poincaré waves (inertia- 
gravity waves). 
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7.3 Gravity waves with centred spatial differencing 


65 


There are several grids known as “Arakawa grids”, which are usually 
identified by the letters A to E (Mesinger and Arakawa, 1976). The 
three most commonly used are illustrated in Figure 7.6. 
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Figure 7.6. The three most common Arakawa grids: A, B and C. 
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For each of these three grids we use the simplest possible centred 
approximations for the spatial derivatives, and when necessary for the 
Coriolis terms. We do not need to study the time differencing since this 


has previously been examined and remains unchanged. 
A-grid: 


Ou; j hi1 j7 hi—ı j 
—— > > as 7.28 
Ot i 2Ax + £3 ( a) 
OVi j hi j+1 h, g—1 
= > > as 7.28b 
Oh; j i+l,j T Ui—1,j ij+1 T Vi, j- 
ail =-H € es Li 4 8 ao 3 +) . (7.28c) 
B-grid: 
ae hag + aa hag — hassi y ftiz, (7.29a) 
x 
Ov; j h; jt T isa j+ 7 hi jT isa j 
— > > > f i 7.29b 
at 4 2Ay Puig ee) 
Oh; j U Ui i—i — Ui-1 Ui—1,j—1 
= H (7 a T 4 7.29 
ðt ( 2Ax ease) 
ae Vig F Uj-1 j ae Uj-1 7) 
C-grid 
Ou; hist — hig 
wid = Th 2 f (Vig + Vitaj + Vitaj- + Vig—1), (7-30a) 
ot Ax 4 
aE = —g E È ! (uij Ui j+ + Uj-1,j41 + Uii), (7.30b) 
oh; ; Ui j — Uiig — Vig — Vig-1 
“H vd Sa vd : 7.30 
ot ( Ag Ay ee 


For simplicity we shall first study the quasi-one-dimensional case 
where u, v and h do not depend on y so that Equations (7.25) reduce to 


Ou T ðh 
at TT Ga’ 
Ov 

OL + fu=0, 
ye a9 


a on 


(7.31a) 
(7.31b) 


(7.316) 
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Inserting the wave solutions from Equation (7.26) into Equations (7.31) 
results in the frequency equation 


() =j oe. (7.32) 


Let us now look at the effect of the finite differencing in space 
for this case. As the variables are assumed not to depend on y, 
Equations (7.28) for the A-grid reduce to 


Ou; ; h; j m h; j 

ae aaa aa a (7.33a) 
Ov; : 
aE — fui j, (7.33b) 
Oh;.; H 

-F = -Ar (Uipi j — Uii), (7.33c) 


and for the B-grid: 


his z pius T Pns hij — Pago pies maa 
Fe = fti (7.34b) 
= se (Uii + Wiji — Uii, — Us-1,3-1) > (7.34c) 

and for the C-grid: 
wes = gas L k / (Vij + Vizag + Vigag—1 + Vij-1), (7.35a) 
Te = = (Uig F Mager + Wiigi Uiig), (7.35b) 
Ta = = (tig = ti-a.) (7.35c) 


Inserting wave solutions with no j-dependence 


(ui, vi, hi) = (Uo, vo, ho el on) 
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Figure 7.7. The function w/f from Equation (7.32), where in the left panel 
gH/(fAx)” = 4, i.e. the Rossby radius is set to two grid lengths 
(VgH/f =2Az) and in the right panel gH/(fAzx)’ = 1/9, i.e. the Rossby 
radius is set to a third of a grid length (\/gH/f = Az/3). The black curve 
corresponds to the analytical solution, the blue curve to the results from the 
A-grid, the red curve to those from the B-grid and the green curve to those 
from the C-grid. NB: The B- and C-grids yield similar results when the 
Rossby radius R is well resolved, but the C-grid results degenerate when the 
grid resolution is coarse. 


(where the imaginary unit is denoted J to distinguish it from the spa- 
tial index 7) into Equations (7.33-7.35) yields the following frequency 
equations for the three grids: 


A grid: (=) =1+ Saen (7.36a) 
9, { Op m gH sin? (kAx/2) 

B grid: ( F ) =1+ P a (7.36b) 
~ (@p\ [kâz _ gH sin? (kAx/2) 

C grid: (2) = cos ( 5 ) a (Ac/2? (7.36c) 


The non-dimensional frequencies wp/f are now seen to depend on the 
two parameters kAx and gH/f?, viz. the Rossby radius R squared, and 
are graphed in Figure 7.7, where they can be validated against results 
from the non-discretised solution of Equation (7.32). 


The pros and cons of the three grids can be summarised as follows: 


e A-grid: The frequency reaches a maximum at kAgx = 77/2, i.e. a 
wave-length of 4 grid intervals. The group velocity is thus zero 
for this wavelength. If inertia-gravity waves of approximately 


The Shallow-Water Equations 69 


this wave number are excited near a point inside the 
computational region, e.g. by non-linear effects or by forcing 
through heating or topography, the wave energy remains near 
that point. Beyond this maximum value, for 7/2 < kAz < 7, 
the frequency decreases as the wave number increases. For 
these waves the group velocity thus has the wrong sign. Finally, 
the two-grid-length wave with kAx = 7 behaves like a pure 
inertial oscillation, and its group velocity is again zero. 

e B-grid: The frequency increases monotonically over the range 
0 < kAx < r. It, however, assumes a local maximum at the 
end of the range, and hence the group velocity is zero for the 
two-grid-length wave with kAx = 7. 

e C-grid: If gH/(fAx)’ > 1/4, the frequency increases 
monotonically in a similar way as in the B-grid case, i.e. when 
the Rossby radius is larger than half a grid length 
(/gH/f > Ax/2). If, however, the Rossby radius is exactly 
half a grid length(y gH if jie 2); the group velocity is zero 
and for smaller Rossby radii the frequency will decrease in an 
unrealistic way with increasing wave number over 0 < kAa < T. 
The advantage of the C-grid lies in that the velocities are 
normal to the grid-box faces, which makes the differencing of 
the continuity equation as well as of the scalar transport in the 
tracer equation a straightforward matter (cf. Section 13.3). 


7.4 The shallow-water equations with leap-frog 


The discretised linearised inviscid shallow-water equations can now be 
written with centred finite differences in both time and space on a 
C-grid as 


n+1 n-1 he = he 
Uij — Mij i+l,j ie oF Eh n i ns 
2At =-g Ar f A (vt, F Use PUT Ufi) , 
(7.37a) 
n+1 =i hv p he 
Vij — Vij ij+1 a S pe i n m 
=-g Uii T Uy ga Ug p41 F Ui 1; 
2At Ay 4 ( ij i,jt+1 t—1,j+1 i aa) , 


(7.37b) 


hett = het yu. — ura f yr = Uaa 
i,j ii H ij ihi y bi i,j ; (7.37c) 
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A von Neumann stability analysis can be applied to the non-rotating 
case, whereby Equations (7.37) become 


À —_ a. elkAx rf 
Shp) 8 ae ho, (7.38a) 
À —_ S ellAy =] 
= h . 
2At Vo g Ay 0; (7 38b) 
À — ATI {= e TkAa 1— e tlAy 


Note that the imaginary unit here is denoted capital “J” in order to 
distinguish it from the index “i”. Equations (7.38a) and (7.38b) can be 
rewritten as 


2gAt (e^ — 1) 


Ug = hem) ho, (7.39a) 
_ 2gAt(e”®^ = 1) 
Up = Ay (A — 3) ho, (7.39b) 


which we insert into Equation (7.38c), resulting in 
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where 


B? = gH(At) eee SA] 


(Aa)? (Ay)? 
Taking the square root of Equation (7.40) results in 


A-A = +41 B, (7.41) 
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which is then transformed into a quadratic equation and has 
the solution 


\ = 421B + V—4B? +1. (7.42) 


The requirement for stability is that |A| < 1, which is satisfied if 
B < 1/2. If we assume that Ar = Ay, then to satisfy the stability 
criterion for all wavelengths we require 


_ VgHAt _ 1 
B= Re Z 4/8 


The dispersion relationship can be found by considering the one- 


zx 0.35. (7.43) 


dimensional case, viz. assuming that the waves propagate along the 
x-axis (l = 0), and substituting 


— „—Iwpât 
A=e CD 


in Equation (7.40) so that 


o 1 . |VgHAt . 
Wp = assin Ae ee) . (7.44) 


The computational phase speed is once again 


Wp 1 A . 
= — = SE ——— 2 A 
Cp A pAg orcsin [2u sin (kKAx/2)], (7.45) 


which is identical to Equation (7.16), which was the computational 
phase speed derived directly from the 1D non-rotating shallow-water 
equations on the staggered grid. 


7.5 Boundary conditions 


There are several types of boundary conditions: Closed boundary con- 
ditions, applied at the border points delimiting land/seafloor from the 
ocean and solid ground from the atmosphere. Open boundary condi- 
tions, taken where the model grid covering the domain under consid- 
eration ends but the real ocean/atmosphere continues. A model can 
also have periodic boundary conditions as previously described for the 
advection equation in Section 4.2. 


7.5.1 Closed boundary conditions 


The staggered B-grid is well adapted to no-slip boundary conditions, 
since the velocity points are located at the corners of the computational 
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grid cell. Unlike the C-grid, there are no ambiguities in the way the 
dynamical boundary condition is imposed at the “corners” of adjacent 
land masses as shown in Figure 7.8. The drawback is, however, that for 
a narrow strait in an ocean model, the B-grid requires at least two grid 
lengths to have a non-zero velocity point as illustrated in Figure 7.8. 
The B-grid yields, however, a more satisfactory dispersion relationship 
than the C-grid since it is better at resolving the Rossby radius at 
coarse resolutions (Batteen and Han, 1981), a feature that makes this 
staggering technique suitable for coarsely-resolved models. 


Figure 7.8. Illustration of how a narrow strait is resolved with a B-grid 


(upper left panel) and with a C-grid (upper right panel). Land cells in 
yellow, ocean cells in white with corresponding u, v and h points in colour. 
The only way to permit velocity points in this B-grid strait would be to “dig 
out” the two cells in green so the strait is at least two grid lengths wide as 
shown in the lower panel. 
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Sponge zone Shallow-water equations 


0 


i Is 
Figure 7.9. Schematic illustration of a sponge zone between the open 
boundary located at i = 0 and i = Ig. 


7.5.2 Open boundary conditions 


An open boundary condition has two main purposes: It should permit 
waves to propagate out from the model domain without being reflected 
back. It should also be possible to force the inner solution with external 
fields, which e.g. can be obtained from observations or models covering 
a larger domain. Open boundary conditions also need to conserve mass 
so that the average of the sea-surface elevation h remains constant. The 
energy budget should also be treated accurately, allowing the correct 
energy flux through the open boundaries to balance the energy flux 
through the sea surface due to the wind stress. 

There are many different types of sophisticated radiative open- 
boundary conditions based on the wave equation. We will here, 
however, only present the simplest, which is the “sponge” boundary 
condition. Here all field variables are first updated using the standard 
interior leap-frog schemes. The field values in the sponge zone are then 
relaxed to the externally given values h? according to 


hH = (1—y) Ay +h", (7.46) 


where h”* is obtained from the model equations, in the present case 
the shallow-water equations. The non-uniformity of h? in the sponge 
zone is taken into account by letting the solution decay as we leave the 
boundary. This decay can e.g. be of an e-folding character or have a 
cosine-shaped relaxation factor such as 


y= 0.5 i + cos (+) (7.47) 


for the interval 0 < i < Is, where J, is the number of grid points in 
the sponge zone, typically 10 to 30. 
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The externally given values h? can originate from observations or 
from another model, which often has a coarser grid. This is the case for 
regional climate models as well as for most local numerical weather- 
prediction models, which are forced at their open boundaries by a 
global circulation model covering the the entire Earth. Figure 7.10 
shows schematically such a nested grid, where the light blue border 
can be treated as an open boundary for the finer interior grid, driven 
by the values obtained from the coarser-grid model The nesting is hence 
a way to “zoom” in on a particular region by increasing the spatial res- 
olution here. The nesting can be one-way, where only the interior values 
are influenced by the exterior values from the coarser surrounding grid 
or two-way, where also the coarser grid values are affected by the data 


from the fine grid. 


t $ 


v 


Figure 7.10. A nested C-grid with a 3:1 ratio between the grid sizes. The 
solid lines denote the coarse grid and the dashed lines represent the fine 
grid. The light blue lines denote the boundary between the fine and coarse 
grids, which can be treated as an open boundary for the fine grid. 
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Radiative boundary conditions, based on the classical Sommerfeld 
condition (Sommerfeld, 1949), were adapted to oceanic and atmo- 
spheric modelling by Orlanski (1976) and followed up by a number of 
other investigators e.g. Higdon (1987) and Flather (1994). These open- 
boundary conditions have the advantage of letting the waves leave the 
domain without reflection and at the same time independently impos- 
ing the open boundary values h”. Tests of different radiative boundary 
conditions can be found in e.g. Nycander and Döös (2003). 


7.6 Conservation of mass, energy and enstrophy 


There are several reasons why numerical schemes for models are often 
formulated so as to respect conservation properties of the governing 
equations. An important practical consideration is that satisfying con- 
servation properties helps to ensure the computational stability of a 
model. Apart from this, the direct physical realism of a conservation 
property may be a desirable feature. For example, ensuring conserva- 
tion of mass prevents the surface pressure from drifting to unrealis- 
tic values in long-term integrations of atmospheric models. Advection 
schemes which satisfy an appropriate dynamical conservation property 
may help to ensure the realism of the simulated energy spectrum. There 
are, however, considerations other than conservation that might influ- 
ence the choice of numerical scheme. Shape-preservation (avoidance of 
the generation of spurious maxima or minima) may be considered as an 
important feature of an advection scheme, and the economy of a method 
(especially the ability to accommodate long time steps) may be a crit- 
ical factor. Indeed, semi-Lagrangian advection schemes (cf. Chapter 
11), originally without formal conservation properties, are increasingly 
being developed for numerical weather prediction. 


7.6.1 The shallow-water equations with non-linear advection terms 


The shallow-water equations with non-linear advection terms will, fol- 
lowing Sadourny (1975), now be presented The momentum equations 
in vector form can in this case be written as 


OV >o r g 
OL +V-VV+fkxV=-gVh, (7.48) 


which can also be expressed as 


ôV > m ls > 
op tE (nV) = (or +5V- 7) , (7.49) 
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where V is the horizontal velocity vector, f the Coriolis parameter, k 
the unit vector normal to the domain S, € = (f + Ov/Ox — Ou/Oy) /h 
the potential vorticity, and h the total water- or air-column height. The 
continuity equation with non-linear terms is 


a +v- (a7) =0. (7.50) 


These equations can also be written in scalar form: 


ðu OB 

oS. (7.51a) 
Ov OB 

a + ¿hu = By (7.51b) 
Oh | O(hu) | O(hv) _ 


where the Bernoulli function is B = gh + Z (u? +v?) =gh ++ 1y VW 
It can be verified that these equations are such that the following prop- 
erties are conserved: 


Total mass: M = i hdS, 
Total Energy: E = f, 4 (gn +7. 7) hdS, 
Absolute potential enstrophy: Z = f 5 5ChdS : 


Here f g dS represents the surface integral over the domain. Enstrophy 
is the integral of the vorticity squared and can be interpreted as a 
quantity directly related to the dissipated kinetic energy. By conserving 
the enstrophy, the model will tend to yield better simulations of eddies. 


7.6.2 Discretisation 


The discretisation on a C-grid is illustrated in Figure 7.11. The spatial 
differencing operators 6, and 6, acting on u, v and h are 


pus Ui,j o 1. § v = of D : (7.52) 
hii; — hi; hisi — hi; 
bh = 1, h= “1G (7.53) 
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Figure 7.11. C-grid with points for the zonal velocity u, meridional velocity 
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v, water- or air-column height h and vorticity ¢. 


Note that all these differences are centred at different points, which is 
best seen from Figure. 7.11. The spatial averages are similarly 


1 1 
u = 3 (tij T Ui-1,;) ; v = z (Viy + Vi j1) ’ (7.54) 
=e. | ee | 
h= 5 (hij + hizig); hk = z (hij + hizt). (7.55) 


It is often simpler to use these differencing and averaging operators 
than to employ index notation. We will use both representations to 
give the reader an opportunity of suffering the agony of choice. 

The mass fluxes U and V are defined at the same points as the 
velocities u and v: 


x 


Ui; =h u = u55 (hig + hi), 
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The gradient operator will act on the Bernoulli function B defined at 
the locations where h is defined: 


2 2 
2 (vz; + ca) 
The potential vorticity is redefined at the corners of the C-grid: 


E= fog ôyu Fig — Vij) (AS = (tiji = uig) /Ay 
vo =7Y g . 


h (hij + hizi + Ragga + Megtg41)/4 


Simple expressions are chosen for a model domain with Nx x Ny 
grid-box faces. 


Nx Ny 
Total mass: M = 90 SUA, Ardy . 


i=1 j=1 


Total energy: E = 47 (gh? + hu?’ + hv?") ArAy 


=i 3 > { ght, + rind [ (aes + u? ay) T (v?, + vj—)] \ Ardy . 


Absolute potential enstrophy: Z = > X. eh” Ardy 
Nx Ny 


ee Ezy alig + Risa + higy + hipiji) ATAY . 


Below the symbol }> refers to a summation of the same species over 
all grid points. Note that due to symmetry 


To = oe 
and that due to skew-symmetry 


X a.b = — X bôza. 


The time derivative of the total energy is 


dE Ou Ov Oh 
ae B . : 
L (v Oe ee ) oO) 
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A simple energy-conserving model can be defined as 


Ou zr 

= +86 B= 
at EV On 0, 
Ov yt 
ee +å, B = 
Ji + EU Oy 0, 
Oh 
Sp +AU +V =0. 


This can also be formulated in a more detailed way using index 
notation: 


a =3 feus (Vig + Viti) + Sij-13 (Vi j-1 + Vanam) 
Bizi; — Bij 
= a 2 (7.58a) 
x 
a = 5 lés (Ui j + Uij+1) + Sizi (Ui—1,; + Use) 
Biju — Bij 
= ey ae a (7.58b) 
y 
Oh; j Ug = Ua Vig = Vai 
= a 7.58 
Ot Ax Ay Mae) 


Energy conservation can be obtained from Equation (7.56): 


Z SE (va - vee) -XC (-U6,B — Bé,U) 
+ X_ (-V6,B — Bô,V) =0, 


where each of the three summations cancel out due to the symmetry 
or skew-symmetry of the operators. 
An absolute-potential-enstrophy model can be defined as 


au zey ih Ha, (7.59a) 
- IFT opi, (7.59b) 
an +6,U +6,V =0. (7.59c) 


In the corresponding vorticity equation, the discretised gradients 
vanish, viz. 6,6, = 6,02, So that 


2 (E) +4. (ET) +4, (@V") =0, 
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which when combined with the averaged continuity equation 
O (= =y" =Y 
SE) +6,(0" ) +a (V) =0 
ai toy 


yields the equation for conservation of the potential enstrophy: 
O =x —ryt =y. 
5 (ER) +6, (EU) +4, (EV) =0. 


7.7 A shallow-water model 


We will here summarise the results above by formulating a model based 
on the discretised shallow-water equations on a C-grid. This will be 
done in the way it is programmed in computer code, and will hence be 
close to how e.g. a Fortran code is structured. The following steps are 
to be taken: 


1. Set the initial condition of the fields for n = 0 over the entire 
model grid indices ¿į and j so that 


n=0 n=0 h=? 
i,j 


Uij s Vij o 


are known. 
2. Integrate the shallow-water equations a first time step with an 
Euler-forward scheme and “loop” over all the model grid indices 


i and 7: 
hPa- h, f 
1 _ 0 i+1,j i,j 0 o 
Uj =u; + At l I— Rr -J (vis + Urea 
Forga + v?) | ) 
h ah, f 
I 240 i,j+1 i,j 0 0 0 
Vij = Uj t At l g Ay A (i: F Ui 544 F Uii jpa 


+u) | , 


0 0 0 0 
üs = Ui Ups Ua 
hi, Z be — AtH ( ij izli y ii ij ) ; 
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Time-integrate the model from n = 1 to n = N,, where N, is 
the total number of time steps to be computed so that the total 
time integration will be N,At. “Leap-frog” the time step with 
loops over 7 and 7: 


h®,,.—h”. f 
n n— i+], 1, n n 
uwy = uns -2At | g oe eg 4 (u F Uti 
Hiti j- t UE 5-1) | , 
hanh, f 
n+1 n—1 ijt 2, n n 
Upp) = Upp + 2At l ge Ay i 1 (up; + ukj 


n n 
FU; a gta + ti 45) | , 


u?.— u”. yr. S 
Ant! = ht — 2AtH ij i—l,j ij ij—1 . 
ij ij Az d Ay 


Apply a Robert-Asselin filter in order to suppress the 
computational mode: 


no- an n—1 n n+l 

Ui j Ugg TY (ur; — 2u; + Uig ) , 
TE pan n-1 n n+l 

Ug = Vag FY (vt; = 20; + Viy ) ) 


hi; = hes + (hiz — 2h; + hey’) f 


Store the resulting fields at regular time intervals and compute 
some statistics, e.g. the total volume V, the kinetic energy Ep, 
and available potential energy Ex: 


Nx Ny 
V=S oS on Ardy, 
i=1 j=l 
oe Ny 
Ep = S DD (M) ArAy, 
i=1 j=1 
H he 


i=1 j=l 
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6. To economise disk space we switch the time-step results (since 
we only store three of these) before returning to the beginning 
of the time loop so that n > n — 1 and n + 1 —> n and 


n-1 _ „n nm—1 _ „n n-1 n 
Ua 5 Ung Vij Vig hij = Mij 
n n+1 n n+1 n n+l 

Ung = Uij o Vij a Mij 5 Baj 


7. End the time loop of the model and the entire model code. 


This shallow-water model can be extended to include terms repre- 
senting non-linear advection and friction/viscosity, the latter to be 
introduced in next chapter. 


8. Diffusion and Friction Terms 


In this chapter we will investigate the discretisation of friction and 
diffusion terms and how this affects the stability of the solution. These 
terms are included in most models from very simple ones based on the 
shallow-water equations to highly complex ocean-atmosphere general 
circulation models. 


8.1 Rayleigh friction 


We start by studying the simplest type of friction parameterisation, 
Rayleigh friction, where the retarding acceleration is directly propor- 
tional to the velocity. A straightforward example is given by 


Ou 
a ery (8.1) 


with the solution 
u (t) = uoe”. (8.2) 


Using a centred time difference (i.e. a leap-frog scheme), which is the 
most common technique employed for equations with advection terms, 
the discretisation of Equation (8.1) is 


n+l n-1 
Uj; 


J = F 
IN = —yuj. (8.3) 


When a stability analysis is undertaken (in analogy with the one the 
advection equation was subjected to in the previous chapter) with 
T= u? A” one finds that 
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If the right-hand side of Equation (8.3) is taken at time step n, 


ie. —yur, then Ai. = —yAt+ 4/1+ (yAt)’, which has at least 
one root that is always greater than one for any yAt > 0. The 


scheme is hence unconditionally unstable. 

If the right-hand side of Equation (8.3) is taken at time step 
n—1, te. —yur then A? = 1 — 2yAt. The scheme is 
conditionally stable since \? < 1 if yAt < 1. But since A? < 0 
for 1/2 < yAt < 1, the roots of A will be purely imaginary and 
the solution u will oscillate and change sign for every second 


time step. If e.g. yAt = 1, then A2 = +i and 

u” = i” =1,0,—1,0,1,... 

For the more restrictive condition yAt < 1/2, A will be real and 
u will have a more realistic evolution in time with no numerical 
oscillations. 

If the right-hand side of Equation (8.3) is taken as an average 
over the time levels n — 1 and n + 1, the following 
finite-difference equation is obtained: 

a a, 


arg te). 


u 


This, as discussed in Chapter 3, is known as the Crank-Nicolson 
scheme and is said to be implicit because it includes a term at 
time level n + 1 on its right-hand side. It yields the best 
approximation of Equation (8.2), and the stability analysis 
results in A? = (1 — yAt)/(1 + yAt) < 1. The scheme is hence 
unconditionally stable. For the same reasons as above one 
requires yAt < 1 in order for a realistic evolution in time. 
Implicit schemes are often complicated to solve since they 
include values on both sides of the equation that need to be 
determined simultaneously. This is, however, not so in this 
particular case, since the right-hand side is evaluated at the 
same spatial grid point 7 as the left-hand side and the equation 
can be rearranged so that 


yrtt = 1—yAt yr 
I 1+~7yAt ? 


When, as in this case, employing the Crank-Nicolson scheme for 
only a time integration, one should use a two-time-step 


o 


Figure 8.1. The Rayleigh friction equation (8.1) integrated analytically 
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(top left) and numerically with the right-hand side of Equation (8.3) at time 


step n (top right), which clearly gives an unstable solution. When integrated 


with the right-hand side at time step n — 1, the solution is stable and smooth 
with yAt = 0.22 (middle left) but oscillating with yAt = 0.8 (middle right). 


Bottom panels show the results of integration with the Crank-Nicolson 


scheme giving stable solutions, but oscillating when yAt = 1.5. 


integration with an Euler-forward scheme so that only two time 


steps are used and the equation becomes 


n+1 __ 


1 — yAt/2 


Ter 
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The stability analysis above is only strictly valid for these discretisa- 
tions of the very simple Rayleigh-friction example given by Equation 
(8.1). However, it turns out that one obtains approximately the same 
stability criterion when a Rayleigh-friction term is included in the 
momentum equations or in the tracer equations in a GCM. It is, how- 
ever, not possible in these cases to undertake a stability analysis of 
these more comprehensive equations. 


8.2 Laplacian friction 


A somewhat more realistic friction parameterisation is based on the 
Laplace operator: 


2 

a ea (8.4) 

Ot Ox? 
where A is the viscosity coefficient with the unit m?/s. The letter A 
originates from the German word Austausch, which means “exchange”, 
referring to the exchange of water “parcels”. It replaces the molecular 
viscosity with a much larger eddy viscosity in order to parameterise 
the sub-grid scales in the momentum equations. Note that often the 
letter K is used in Equation (8.4), known as the heat equation, when 
it represents the diffusion of a tracer. This equation is known as the 
diffusion equation and is a parabolic PDE. 

For a single wave number k Equation (8.4) has the solution 


u(x, t) = uyet*e- Ae", (8.5) 


The simplest way to construct a finite-difference approximation of a 
second-order derivative is to apply finite differencing to a finite differ- 
ence. This is achieved by first postulating two finite differences centred 
on the intermediate positions j + 1/2 and j — 1/2 as illustrated by 


Figure 8.2: 
du Uj — Uj—-1 
— x — 8.6 
(=) _ Ar 7 PE 


=) Uj41 T Uj 
ae ga (8.7) 
(= peta Az 


Since the second-order derivative is defined as the derivative of the 


? 


? 


derivative, we can similarly construct a further finite difference: 
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TA 

du _ d (du ea = aoe _ Wy 2u; + uji 

dx? j ~ dx \ dx Ax (Az)? 


Figure 8.2. The second derivative is the derivative of the derivative. By 
first estimating the finite differences at j + 1/2 and j — 1/2 and then the 
finite difference of those two, one obtains the second finite difference at 7. 
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Figure 8.3. The heat equation integrated numerically using Equation (8.12) 
with the right-hand side at time step n and with v = 0.01 (left panel) and in 
the right panel with the right-hand side at time step at n — 1 and v = 0.125. 


eu... aes R oem ee Oua Fu 
dr? jJ, |dx \dz/], Ag (Ax)? 


(8.8) 
The advantage of this formulation is that it is straightforward and 
intuitively evident. The disadvantage is that it does not provide an 
estimate of the accuracy of the scheme. To obtain this we use the 
Taylor-series method previously employed in Section 3.2. A centred 
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finite difference of the Laplace operator corresponding to the second- 
order derivative can hence be obtained by combining two Taylor series: 


7 du (Az)? (du (Az)? (du 
a ad (a) + 2 (= F = 6 dx? } , 


8.9 
24 dx‘ }, T? 
a A du\ | (Az) (Pu (Az)? (@u 
S (æ), 2 dz? j 6 dz? j ae 
j T i (8.10) 
(Ax) (dtu 
"24 dzt j ae 
3 
Adding these two equations and dividing by (Ax)? we obtain 
Uj+ı 7 2u; + Uj-1 du 1 2 dtu 
= A = sits 8.11 
(Ax)? (oa j i 12 yaa dzt } , + oaa 


This finite-difference approximation of the second-order derivative is 
hence accurate to order (Az)’, which is the same as saying it has a 
second-order truncation error. 
The analytical heat equation (8.4) can now be approximated by 
integrating in time with a leap-frog scheme: 
i a quit T 2U; + Uji 
2At (Ax)? 


(8.12) 


or 
upt! = a" + 2y (544 EN 2u; + U;—1) ) (8.13) 


where v = AAt/ (Ax)? is the non-dimensional von Neumann number, 
sometimes also called the diffusion number. The von Neumann number 
is, as we will see, now a number that should sometimes not be exceeded 
in order to have numerical stability when integrating an equation with 
Laplacian diffusion. 

A stability analysis using the von Neumann method is undertaken 
by inserting u? = uA” eise 
results are obtained depending on at which time step the right-hand 
side is chosen. Let us examine the same three cases as we did in the 


into this equation. Different numerical 


Rayleigh-friction example above: 
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If the right-hand side of Equation (8.12) is taken at time step n, 
the equation for the amplification factor becomes 


kA 
A? + 8v sin? (=) A-—1=0, 


which has the roots A,» = —a + Va? + 1, where 

a = 4v sin? (422), For the second root it is recognised that 

A2 < —1 for any v > 0, implying that the scheme is 
unconditionally unstable. 

If the right-hand side of Equation (8.12) is taken at time step 
n — 1, the amplification-factor equation becomes 


A 
X = 1 — 8v sin? (=). 


The scheme is stable when —1 < A? < 1, which is the case when 
v < 1/4, and thus the scheme is conditionally stable. However, 
for the same reasons as for the Rayleigh-friction equation, we 
recommend the stricter condition v < 1/8, this in order to have 
à? > 0 and hereby avoiding oscillations in time of the solution. 
If the right-hand side of Equation (8.12) is taken as an average 
of time levels n — 1 and n+ 1 (the Crank-Nicolson scheme) we 
have 


ust! — up! _ A (usti — 2upt* + upt | ure — Que + ue 
2At a2 (Az) (Ar) 


(8.14) 


This scheme is implicit as it includes terms at time step n + 1 
on the right-hand side of the equation, which, however, can not 
be solved as easily as in the Rayleigh-friction case, this since 
the n + 1 terms on the right-hand side occur at the spatial grid 
points j — 1, j, j + 1. It is, however, possible to use Gaussian 
elimination to deal with these terms. We can nevertheless 
undertake a stability analysis and calculate the amplification 
factor, which is found to be 


2 _ 1—4vsin® (kAa/2) 


= i 8.15 
1+ 4v sin? (kAx/2) ae 
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Here the right-hand side is always smaller than one and the 
scheme is hence uncondionally stable. In order to avoid 
imaginary roots that lead to oscillating solutions one should, 
however, use v < 1/4. 


In most cases when modelling the atmosphere or the ocean, y and A are 
of such magnitudes that the stability criterion derived in the present 
chapter permits At to be much larger than the value conforming to 
the CFL criterion, which requires a restriction of the Courant number 
u = cAx/At. A common mistake when writing a simple model code is, 
however, to use the unconditionally unstable scheme with the friction 
taken at time step n. 

Note that the schemes in the two last cases discussed above are in 
fact two-level schemes, since we do not use any values at time step n, 
but only at n — 1 and n+ 1. There is consequently no reason to use a 
leap-frog scheme here, and we can instead use an Euler-forward scheme 
in time and replace all time levels n — 1 by n. The stability analysis 
remains unaltered, but, since the time step is halved, we should replace 
At by At/2. It is nevertheless easier to demonstrate the differences 
between the three cases by using leap-frog schemes for all of them. 

In Section 4.3 we saw that wave propagation with the discre- 
tised advection equation required restrictions on the Courant num- 
ber u. Here, we have seen that similar stability conditions arise when 
Rayleigh and Laplacian friction are used. This leads to restrictions 
on the non-dimensional number yAt and the von Neumann number 
v = AAt/(Az)’. The discretised momentum equation will need to 
satisfy all these stability criteria when friction parameterisations are 
included and the CFL criterion is satisfied. The next section will exam- 
ine how the stability criteria associated with advection and diffusion 
are interrelated. 


8.3 The advection-diffusion equation 


Let us now examine an equation with both advection and diffusion 
terms: 
Ou Ou aoe 
Bt On ~ Oa?” 
which in this form combines both the parabolic and hyperbolic prop- 


(8.16) 


erties of a partial differential equation. It has the analytical solution 


u (x, t) = uge EE-D At, (8.17) 
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We have previously seen that a discretisation of the advection equation 
with a scheme centred in time as well as in space is stable, while for the 
diffusion equation the discretised Laplace operator must be taken at 
time step n — 1 in order to ensure stability. Let us now combine these 


schemes: 
ue u Rate -uja L vith — 2u; + uni} (8.18) 
2At i 2Ar (Az)? l 
or 
ut) — uP? — (ut, wa) + 2v (ur — 2u? uti). (8.19) 
We now undertake a stability analysis by inserting u? = A” e*IA® into 
this equation, which after some calculations yields 
A? + iar +8b—-1=0, (8.20) 


with the coefficients a = usin (kA) and b= v sin? (kAw/2). 
The solution of Equation (8.20) is 


A= —ia + V—-a? + 1 — 8b. 


If 1 — 8b — a? > 0 then |A|? = a? +1 — 8b — a? = 1 — 8b < 1, viz. the 
scheme is stable for this root. 

If 1 — 8b — a? < 0 then \ = —i (a F Va +8b-1) > X=- 
(2a? + 8b — 1 F 2av/a? + 8b — 1). 


It is not immediately evident when the second root of this expression 


for à? yields a stable solution, and thus we have graphed \? as a function 
of a and b in Figure 8.4. 

The stability analysis above can be tested by reformulating Equation 
(8.19) as 


n+l 


Uj 


=u — u (u w a) + 2v (uy — 2u un), (8.21) 
which we then integrate numerically for 100 time steps with the same 
initial condition as for the Rayleigh and diffusion equations in the pre- 
vious sections. The choice of both the Courant and von Neumann num- 
bers will determine the stability of the integration. From Figure 8.4 we 
can see that a stable and non-oscillating solution will require a Courant 
number u below 1. If we choose e.g. v = 1/16, we recognise from Figure 
8.4 that |A|? < 1 for u up to approximately 0.70. To test this we have 


integrated Equation (8.21) with a Courant number u just above and 


92 


Basic Numerical Methods in Meteorology and Oceanography 


Unstable 


0, 
Sy 
Oy 
Lo. 


%; 
o 


00 01 02 03 04 05 0.6 07 08 09 10 11 


-1.2 -1.0 —.80 —.60 —.40 -.20 0 .2(0 .400 .600 .600 1.00 1.20 


Figure 8.4. The squared amplification factor À? as a function of the 


Courant number u and the von Neumann number v. The purple region, 


where \? < —1, corresponds to where the solutions are unstable. The blue 


region, where —1 < Al? < 0, represents stable but oscillating solutions. The 


red region is for 0 < A? < 1, which corresponds to stable solutions with no 


oscillations. The two yellow dots indicate v = 1/16 with u = 0.70 and 


u = 0.76, which are the two test cases illustrated in Figure 8.5. 


below this critical value. Figure 8.5 shows the results of these two inte- 
erations, with a stable solution obtained for u = 0.70 and an unstable 
one for u = 0.76. 


Exercises 

1. Undertake a stability analysis for the Rayleigh-friction equation 
with the right-hand side of Equation (8.3) taken at time step n. 

2. Same as in 1) but for the right-hand side taken at time step 
n-li. 

3. Same as in 1) but for the right-hand side taken at time step 
mad: 

4. Calculate the stability criterion for 


ðu o?u 
ashe eg 
fC ep 
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Figure 8.5. The heat-diffusion equation integrated numerically using 
Equation (8.21) with v = 1/16 and u = 0.70 in the left panel and u = 0.76 
in the right panel. Note the visibly growing instabilities in the right panel 
after 60 time steps. 


using the following scheme: 


yee -u7 — juin — 2up + uly 
2At (Ar)? 


Estimate an upper limit for At when 
i) A = 10°m?/s, Ax = 400 km (large-scale horizontal 
diffusion), 
ii) A = 1m?/s, Ax = 10m (vertical diffusion in a boundary 
layer). 
5. The diffusion equation can be integrated using the 
Crank-Nicolson scheme: 


THT AJI a TT, r nrn aT 


A 2 (Ax)? (Ax)? 


Examine the stability of this scheme! 


9. The Poisson and Laplace Equations 


Consider the elliptic Poisson equation in two dimensions: 


o? o? 
V’u = (= + Z) u = f (x,y). (9.1) 


If f(x,y) = 0 this is known as the Laplace equation. As shown by 
Figure 9.1, Equation (9.1) can be discretised on a grid: 


Ui-1,j — Ui j F Uitij | Uig- T Ui j H Uj 


(Ax)’ (Ay)? 


= Jij) (9.2) 


which can also be written as 


(Ay)? (Uii, + Uist) + (Ax)? (Ui j-i + uijt) — (ArAy)” Joz 


k 2 [(Aw)* + (Ay)”] 
(9.3) 
If we consider a square grid such that Ax = Ay, Equation (9.3) sim- 
plifies to 
1 2 
t= [iij + Uiga + Uig- Fijt — (Az)” fagl: (9.4) 


When the boundary values for the domain are known, it is possible to 
solve this finite-difference equation by iteration. 

In iterative methods we need initial values at iteration level uj", 
m+1 
ij 
is repeated until the difference between the results from two successive 


(m = 0 initially) and the purpose is to calculate u; . This procedure 


iterations decrease below a prescribed value at each grid point. 
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1 i-] i i+] IX 


Figure 9.1. Grid for the Poisson and Laplace equations. Boundary values 
required for the four walls at the blue points. The green points illustrate 
which points are needed to compute the red point. The solution is obtained 
by calculating iteratively starting from the orange arrow. 


9.1 Jacobi iteration 


Values from the previous iteration level are used, which results in 


up = Z ane F Uiga, gj T Uig- F Ugg T (Az) fasl- (9.5) 
The method works but is somewhat inefficient and is not used for solv- 
ing practical problems. 
A Fortran-code segment of this iterative computation could look like 
this: 


u=0. ! initialise the field to zero 

f=0. ! Set the function f to something 
omega=1.5 ! relaxation factor 

! Set the boundary conditions to something 
u(i,:,:)=10. ; u(IX,:,:)=10. 

u(:,1,:)=10. ; u€:,JY,:)=10. 

do m=1,500 ! number of iterations 

do j=2,JY-1 


! Jacobi iteration 

! u(i,j,m+1)=0.25*(u(i-1,j,m)+u(i+1,j,m) & 

! & tuli jism) tulis jt+1,m)-dx**2*f(i,j)) 
! Gauss-Seidel iteration 

! u(i,j,m+1)=0.25*(u(i-1,j,mti)+u(iti,j,m) & 

! & tuci,j-1,m+i)+ucdi,j+1,m)-dx**2*f(i,j) ) 
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! Successive Over Relaxation (SOR) 
u€i,j,m+1)=(1.-omega)*u(i,j,m)+ & 
& omega*0.25*(u(i-1,j,m+i)+uCitl,j,m) + & 
& u(i, j-i,mti)+u(i,j+1,m)-dx**2*f(i,j) ) 


res(m)=res(m)+(uC(i,j,m+1)-u(i,j,m)) **2 
enddo 

enddo 

print *,m,res(m) 

enddo 


9.2 Gauss-Seidel iteration 


A clear improvement in efficiency of iterative methods is obtained if 
we use the newly computed values in the iteration formula: iteration 
level m + 1 values are available for nodes (i — 1,7) and (i, j — 1) when 
calculating u for node (i, j). Thus the Gauss-Seidel formula is: 


ti = a Par F Unga + u + Ui jp T (Az) fasl: (9.6) 
The inclusion of the two newly computed values makes Gauss-Seidel 
iteration more efficient than Jacobi iteration. Note that the Fortran 
code above for the Jacobi method can easily be “upgraded” to the 
Gauss-Seidel method by simply removing the iteration index m. 


9.3 Successive Over Relaxation (SOR) 


The Gauss-Seidel iteration method can be further improved by 
increasing the convergence rate using the method of Successive Over 
Relaxation (SOR). The change between two successive Gauss-Seidel 
iterations is denoted the residual c, which is defined as 

C=a E T (9.7) 
In the method of SOR, the Gauss-Seidel residual is multiplied by a 
relaxation factor w and a new iteration value is obtained from 


=u +we = ur +o (ar —um) = (w) wie, (9.8) 


m+1 
j ij i,j 


Uig 


where oe denotes the new iteration value obtained from the 


Gauss-Seidel method using Equation (9.6). It can easily be seen 
that if w = 1, SOR reduces to the Gauss-Seidel iteration method. 
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Figure 9.2. Numerical solution of the Laplace equation using the supplied 
Fortran code. Iterated 500 times using the Jacobi method (upper left), the 
Gauss-Seidel method (upper right), the SOR method with w = 1.5 (lower 
left) and with w = 2 (lower right). Note that this last case does not converge 
because of too large a relaxation factor. The end solution should converge 
towards 10, since that is the value of all the boundary conditions and f = 0. 


By substituting Equation (9.6) of the Gauss-Seidel iteration method in 
Equation (9.8), we obtain the equation used in the SOR method: 


m m wW m m m m 
u =(1- w) uij t a [a + Uiga, F up +U T (Ac)? fasl: 
(9.9) 
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Usually the numerical value of the relaxation factor w can be obtained 
by trial and error and the optimum value is generally around 1.5. In 
the case that 0 < w < 1 , the method is said to be “under-relaxed” . 
According to the choice of the parameter w, we either extrapolate for 
w > 1 or for 0 < w < 1 interpolate between the old iteration value at 
level m and the Gauss-Seidel value at level m + 1. If we extrapolate 
too much, i.e. w is taken too large, the iteration starts to oscillate and 
probably collapses. The iterative methods described above are often 
referred to as relaxation methods as an initially guessed solution is 
allowed to slowly relax, reducing the errors, towards the true solution. 

Finally, it also deserves mention that multi-grid methods can also 
be employed, where sequences of coarser grids are used so as to provide 
initial guesses for the finer grids. This is done in order to speed up the 
convergence of the iterative procedure. See e.g. Hackbusch (1985) for 
a comprehensive overview of multi-grid methods and applications. 


9.4 Helmholtz Decomposition 


An example of Poisson equations, where one uses the iterative meth- 
ods above in order to solve the equation, is when computing stream 
functions and velocity potentials. The Helmholtz theorem states that 
a velocity field can be decomposed into rotational and divergent parts: 


U = Uy + Uy, (9.10a) 
v = Uy + Vy, (9.10b) 


where the subscript x denotes the divergent irrotational part and w the 
non-divergent rotational part. The stream function is defined as 


Op Oy 


Ae = Vy By = —Uy- (9.11) 


The velocity potential is defined as 


ôx ôx 


De =the ge =e (9.12) 
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An equation for the stream function Y% can be derived by 


09.10b 9.10a_ 

Ox Oy ` 
ðv ðu Ovy | ? _ 
ðr Oy Ox Ox Oy Oy (aaa 
Ow ; 0?x oy OX o 


where € is the relative vorticity. 
An equation for the velocity potential can similarly be obtained by 


Or Oy (9.14) 
Vp=V-V. 


9.10a  ə9.10b 


These two Poisson equations may be solved iteratively as explained 
previously in this chapter. 


Exercise: 


Set up a numerical model of the Laplace equation with a grid of 10 x 
10 points. Start with u = 0 in the interior and with u = 1 as boundary 
conditions. Test the convergence of the 3 different iteration schemes 
from this chapter. 


10. Implicit and Semi-Implicit Schemes 


The time step permitted by the economical explicit schemes, twice 
that satisfying the CFL criterion, is still considerably shorter than that 
required for accurate integration of the quasi-geostrophic equations, 
which do not permit fast oscillating waves. Thus we will here consider 
implicit schemes, which have the pleasing property of being stable for 
any choice of time step. 


10.1 Implicit versus explicit schemes, a simple example 


For implicit schemes, the spatial terms are evaluated, at least partially, 
at the unknown time level. Let us consider one of the simplest possible 
examples by examining the one-dimensional diffusion equation, also 
known as the heat equation (8.4): 

Ou Oru 

eae (10.1) 
A formal solution of this parabolic differential equation requires an 
initial condition as well as two boundary conditions, the latter, in order 
not to complicate our problem unnecessarily, here taken to be Dirichlet 
conditions (cf. Section 2.1) 

This equation is discretised with centred spatial finite differences 

and integrated in time with an Euler-forward scheme. In traditional 
explicit form we obtain 


2u? + uka), (10.2) 
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which can be rewritten as 


up’ = vuz + (1 — 2v)u; + Yury (10.3) 


where v = AAt/ (Ax)? is as previously the von Neumann number. This 
equation is explicit in terms of u”*', which is the value at the unknown 
time level n+ 1, and is hence possible to solve. A stability analysis can 
be performed and shows that it is conditionally stable (—1 < A < 1) 
for v < 1/2. A stricter condition with only a positive root (0 < A < 1) 
for the non oscillating solution is obtained for v < 1/4. 

A similar approach, but evaluating the spatial term at the unknown 
time level n + 1, yields the fully implicit discretisation 


upt = up + njutt — duet? + unt), (10.4) 


which can be rewritten with all terms at time step n+1 on the left-hand 
side as 


— vwt + (1+ 2v)uph — vu = up. (10.5) 
This implicit discretisation is unconditionally stable. To solve the 
equation one needs to consider all grid points 7. In the present case, 
when we are dealing with the linearised heat equation, the problem can 
be expressed as a linear system of equations AX = B, where A isa 
matrix, X a vector given by the unknown values of u at time n + 1, 
and B a vector given by the known values of u: 


(1 + 2v) =v uz" uz + vut" 
=v (1+ 2v) =y uzt! uz 
—v (1+ 2v) -y unt] ul 
—y (1+2v)| (upt un_, teu 
(10.6) 


For didactic reasons we take u7*' and u?** to be known from Dirich- 


let boundary conditions. Neumann and Cauchy conditions can equally 
well be applied, but are somewhat more complicated to implement. 
The solution at time level n + 1 is determined by solving this system 
of equations. The implicit method is consequently very computationally 
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demanding compared to the explicit method, but since it is 
unconditionally stable it is possible to use larger time steps. In the 
present case, the matrix is tridiagonal, which is advantageous from a 
computational standpoint, since the problem can be solved using e.g. 
the Thomas algorithm, a simplified version of Gaussian elimination. 


10.2 Semi-implicit schemes 


Semi-implicit schemes evaluate the spatial derivative at an average of 
the time levels n and n + 1 instead of only at n+ 1 as in the fully- 
implicit case. If F'(x,y,t) is a term comprising spatial derivatives of 
a given scalar T(x, y,t), we can consider the general expression for a 
discretised version of the equation for the time evolution of u; 3: 


du ul —U", 
=F = 2) bd (1 F”, + BFM 10. 
dt (x, y) At ( p) 1,9 p 1,9 ’ ( 0 7) 


where 8 = 0 yields an explicit scheme, 8 = 1 a fully implicit scheme 
and 0 < 8 < 1a semi-implicit scheme. 

A commonly used semi-implicit method is given by the Crank- 
Nicolson scheme, in which 6 = 0.5 and the time derivative is expressed 
with the usual Euler-forward scheme. The term comprising spatial 
derivatives is therefore centred at time level n + 1/2, which in fact 
turns this scheme into a trapezoidal implicit scheme in time. By car- 
rying out a Taylor expansion around (i,n + 1/2), one can verify that 
this scheme is characterised by a second-order accuracy in time, which 
represents an appreciable improvement with regard to the first-order 
accuracy of the Euler-forward explicit scheme. 


10.2.1 The one-dimensional (1D) diffusion equation 
The heat equation (8.4) is usually associated with centred differencing 
in space. When using a semi-implicit time scheme it becomes 


n+1 n n n n+l n+l n+1 
i. UR Ui, — 2U; + uy — 2u; + uii 


4 i _ Uit1 
a = BA + ga y 


(10.8) 
The semi-implicit Crank-Nicolson scheme (8 = 0.5) results in a numer- 
ical precision of second order both in time and space and hence the 
truncation error is of O [(At)? ; (Aa)*|; 
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10.2.2 Two-dimensional (2D) pure gravity waves 

Let us now discretise the equations for two-dimensional shallow-water 
gravity waves, viz. Equations (7.25) without the Coriolis terms, using 
the Crank-Nicolson scheme on a C-grid and an Euler-forward time 
integration: 


nmr n gAt n n n n 

wj =u 2AT (hiiia — hag + hirig — his) (10.9a) 
n+l _ „n gAt n+l n+1 n n 

Vij S Uii T 2Ay (Mori = hij thijm hes) (10.9b) 
n+l n 

he = pr, (10.9¢) 

HAt at wti FUZi — UF ay | vp’ vpt H Ui; = Vij- 

2AT 2Ay 


(10.9d) 


A stability analysis of these equations can be undertaken by inserting 
the wave solutions 


(u”, v”, h”) = (uo, tal re he, (10.10) 
and we find 

At 

ig a= a An Cay) (ee ay, (10.11a) 
A 

Up (1 — A) = T (1+ A) (e”*¥ — 1) ho, (10.11b) 

J= e TkAx t= e TlAy 
1-A)=HA . . 


By eliminating uo, vo and ho, we obtain the following quadratic 
equation for A: 


i= 
Mao E 6: 
ra (10.12) 
in? (kAx/2) sin? (lAy/2) ` 
pennar 
oa l Ax? ü Ay? l 
with the two roots 
1-B+2VB 
Va (10.13) 


1+B 
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These two amplification factors have the absolute value 


> (1=-BY VB \ 

|Ax2| = (=) + 14B =]; 
and thus the scheme is unconditionally stable. The scheme is also said 
to be “neutrally stable” since |A| = 1 is just at the edge of stability. 
This example of an application of the Crank-Nicolson scheme shows 
the power of semi-implicit methods as these both decrease the tem- 
poral truncation error from O(At) to O [(At)"] and make the scheme 
unconditionally stable. 

To solve the system constituted by Equations (10.9), as will be done 
in next section, the quantities 6,u”*! and 6,v"*! can be eliminated 
from the third of these equations by applying the operators 6, and 6, 
to the first and second, respectively, and substituting the results into 
the third equation. This yields a set of equations for the height h in the 
form of a linear matrix system involving each grid point of the domain. 
This problem can be solved using e.g. an iterative procedure similar to 
those further discussed in the previous chapter: 


1. Make a first guess h”*? which is usually h”. 
At each of the grid points the value of h”*! has to satisfy the 
equation. 

3. The preceding step is repeated as many times as needed to 
make the change at every point less than some pre-assigned 
small value. 


10.3 The semi-implicit method of Kwizak and Robert 


When considering the shallow water equations Kwizak and Robert 
(1971) chose to use the leap-frog scheme and a semi-implicit differ- 
ence system for variables at time level n + 1. The governing equations 
can be written in a compact form: 


Ou Oh 


ye 10.14 
or Iz Au, (10.14a) 
Ov Oh 

meee Lan 10.14 
Al 5, Av, (10.14b) 


L, (> + m) F An, (10.14c) 
Ot y 
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where A,, A, and A, represent terms that were omitted in Equations 
(7.25) describing the propagation of pure gravity waves. This time we 
apply implicit differencing over a time interval 2At centred around time 
n for the terms containing spatial derivatives by using 8 = 0.5 with 
time steps n—1 and n+1 rather than n and n+1 as previously. Second- 
order centred schemes are used for spatial derivatives and the leap-frog 
scheme for the time derivative, and hence the discretised system is 


urt! = yr! = gAt (oan 4 6h") + 2AtA”, (10.15a) 

ut =u” — gAt Own + Og") +2AtA”, (10.15b) 

Arti = het — HAt (6,u"-? Ẹ ĝu”! ae 6,untt 4 ô v” Tt) + 2AtA}. 
(10.15c) 


We now apply the operator 6, to the first and 6, to the second of these 
equations, and add the results. By introducing the notation 


Out = ôs (sh) and ph = ô, (ô h), (10.16) 
we obtain 


(6,u + 6,v)"** =(6,u + ôv) 
= gAt [See + dyy) A + (See + Syy) AT] 
+ 2At (ô Au + 6,A,)”. 


Substituting the right-hand side into Equation (10.15c), and defining 
the finite-difference Laplacian by 


V2 = bee + Syy; (10.17) 
we find that 


Ant) = A"! — QHAt (6,u+ ôv) T + HAP (VZR! + Vh) 
+ 2At[A, — HAt(5,A, + 6,A,)]”. (10.18) 


By, in addition, introducing the definitions 
F° = ph"! — 2HAt(6,u+6,v)” + gHACPVZA", — (10.19) 
G” = 2At [A, — HAt (ôs Au + 6,A,)]”, (10.20) 
this result can be formulated as 


net! — gHAPV?h H = F” 4G”, (10.21) 
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where the terms have been arranged to show that at time level n, the 
right-hand side of this equation is known at all grid points. This is an 
elliptic PDE reminiscent of the Helmholtz equation 


VW h+ah+b(z,y) =0. (10.22) 


Several methods are available for resolving this standard problem. Once 
it has been solved for h”+!, then u”+! and v"t! can be obtained directly 
from Equations (10.15). 


11. The Semi-Lagrangian Technique 


In an Eulerian advection scheme an observer at a fixed point in point in 
space watches the surroundings. Such schemes work well on structured 
grids of the type to be discussed in Chapter 12, but often lead to unnec- 
essarily restrictive time steps imposed by the requirement of computa- 
tional stability. In a Lagrangian advection scheme the observer watches 
the ambient world evolve while travelling “on board” a fluid particle. 
An advantage with a Lagrangian scheme is that one can use much 
larger time steps than in the Eulerian case. A disadvantage is, however, 
that a regularly spaced set of particles will in most cases subsequently 
evolve into one which is highly irregularly spaced and important char- 
acteristics of the flow may consequently be lost. The advantage with 
the semi-Lagrangian advection schemes is that they combine the reg- 
ular resolution of the Eulerian schemes with the enhanced stability of 
the Lagrangian ones. Robert (1981) proposed using a semi-Lagrangian 
scheme for the treatment of the advective part of the equations (see 
e.g. Staniforth and Côté (1991) for a general review). 


11.1 The 1D linear advection equation 


To present the basic idea underlying the semi-Lagrangian method in 
its simplest context let us first examine the one-dimensional advection 
equation formulated within an Eulerian framework: 


OF OF 
ai, wa 

ðt Be ou 

where c is a given constant velocity and F the passively advected tracer 


such as e.g. the temperature or moisture. The most straightforward 
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-4 j-3 j-2 j-1 j jtl jt2 j+3 j+4 


Figure 11.1. The green stable region for the advection equation discretised 
with centred finite differences in both time and space. The maximum 
distance you may travel during one time step At is one grid length Ax. The 
red dotted line illustrates how advection with a speed c = 4A2/3At takes 
place in the red unstable region, which can, however, be solved numerically 
with a semi-Lagrangian scheme. 


discretisation of this equation, as presented in Chapter 4, is centred 
differencing in both time and space: 
n+l n—1 n n n 
o eae (11.2) 
2At 2Anz 

This scheme will only yield a stable solution when integrated with a 
Courant number u < 1 as shown by the green area in Figure 11.1. 

Equation (11.1) can be formulated in a Lagrangian framework 
instead of an Eulerian one, resulting in 


DF _ 
Io 
where the Lagrangian derivative is defined as D/Dt = 0/0t + cO/Oz. 
Equation (11.3) simply shows how the value of F is constant along the 
corresponding trajectory. 
In discretised space this implies that pot must be equal to the value 
of F at time step n, which can be expressed as 


0, (11.3) 


Pett = pe (11.4) 


where * symbolises the spatial position at time level n, which is nor- 
mally not a grid point and in our case is where the red dotted line 


The Semi-Lagrangian Technique 111 


in Figure 11.1 crosses time level n (shown as the red circle located 
between j — 2 and j — 1). The value of F? can hence be obtained by 
interpolation between these grid points: 


Fe =oF",+(1—a)F",. (11.5) 


Here a = frac (u), where 0 < a < 1, is the fractional part of the 
Courant number. Making use of the integer part of the Courant num- 
ber, p = int (u), viz. u = p + a, this relationship can be expressed in a 
more general way as 


Ft =aF", ,+(1-a)F 


—p* 


(11.6) 


In our example p = 1 and a = 1/3, as shown in blue in Figure 11.1. The 
discrete expression can be obtained from Equations (11.4) and (11.6), 
resulting in 

Fe} —@F" 4. +(1—o0) F",. (11.7) 
Figure 11.2 shows how this semi-Lagrangian discretisation of the linear 
1D advection equation can be used with = 4/3 (left panel) but is 


clearly unstable (right panel) for the centred scheme in time and space, 
which was presented in Chapter 4 and Figure 4.4. 


uunulonnnshaiib a tnabar banteikarnren 
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J 
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J 
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Figure 11.2. Hovmöller diagrams of the linear advection equation 
integrated numerically with u = 4/3. A semi-Lagrangian scheme has been 
used for the left panel and a centred scheme in both time and space for the 
right panel. Note that the integration with the centred scheme blows up 
after just a few time steps, cf. Figure 4.4. 
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11.2 Stability analysis 


A von Neumann stability analysis is undertaken of the semi-Lagrangian 
discretisation of the 1D linear advection equation discussed above. We 
search for a solution of the form F? = Ar Foe**J4* that we substitute 
in Equation (11.7), which results in the amplification factor 


d= ace OTDA + (1 — ae PA? = erse (1 — a + ae ****) . (11.8) 
The scheme is stable for |A| < 1, which is why we consider 
AP = ie rae (1 —atae‘**) |’ ; (11.9) 


We use Euler’s formula, et? = cos 8+isin 6 and the associated absolute- 
value result |cos 6 + isin 8| = 1, which leads to 


||? = 1 — 2a(1 — a) [1 — cos(kAz)]. (11.10) 


The minimum value of this expression is obtained when cos(kAxr) = 
—1, which yields 


|A|? = (1 — 2a)? < 1. (11.11) 


The maximum is obtained for cos(kAz) = 1: 
A? =1. (11.12) 


The scheme is hence unconditionally stable. The time step can thus be 
much larger than in the case of e.g. the leap-frog scheme. A larger time 
step will, however, decrease the numerical accuracy in a GCM, a fact 
that must be taken into account. In a GCM employed for numerical 
weather prediction this will nevertheless make it possible to use a time 
step approximately six times larger than in the leap-frog case. The 
absolute stability of the semi-Lagrangian scheme can be understood in 
the sense that taking one single step along the flow in both time and 
space is a way to adjust the spatial resolution. This is like prescribing 
the spatial resolution coarser when possible, “as if the Courant number 
had been equal to one”. 


11.3 The advection equation with variable velocity 


The derivations in the previous section can be extended to the 1D non- 
linear advection equation for a tracer, where we replace the constant 
phase speed c with a velocity, which varies in time and space. 
DF OF | OF 
= u 
Dt ðt Ox 


=0, (11.13) 
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where u = u(x,t) is a given velocity and F the passively advected 
tracer such as e.g. the temperature or moisture. It can also, as in the 
momentum equations, be the velocity itself so that F(x,t) = u(x,t) 
We consider the case where the velocity u is known at all grid points 
in space and time. The centred scheme in both time and space in the 
Eulerian framework is then 

Ft — ae Fe 


| u”? mr ha 
2At © Í INT 


=0. (11.14) 


The Courant number is now u = u% At/Ax and hence variable in time 
and space but should never exceed 1 anywhere and anytime for this 
centred case. Figure 11.3 shows an example of this, where the velocity 
has been prescribed as varying in both time and space, and where the 
maximum velocity was just at the limit of stability (u = 1). 

The semi-Lagrangian scheme for a variable prescribed velocity is, 
just like in the linear case, based on that one can compute the value 
F?+! by following a trajectory backwards from its upwind position at 
a previous time step. The “exact” backward trajectory from the point 
where F’"*? is illustrated by the blue curve AC in Figure 11.4. For 
our calculations it will be approximated by the red straight line A’C. 
Equation (11.13) states that the scalar F remains constant along a fluid 
path or trajectory. The integration along the approximated trajectory 
of Equation (11.13) is thus 


PeH = Fo (11.15) 


where F'”~' is located at A’, which subsequently is to be determined. 
The particle displacement in the x-direction over the two time steps 
from point A’ to C is 2u”At. Here u” is the interpolated velocity at 
point B: 


us = au; p1 l= ays (11.16) 


where p = int(u) and a = p — p. This is the same interpolation as 
for the case with a constant velocity but with the difference that the 
Courant number now depends on a variable velocity u”: 


_ u%At 
a Ag ` 


(11.17) 


Since both a and u? depend on each other we need to iterate in order 
find the point B with the corresponding velocity wu”. Once this is found 
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Figure 11.3. Hovmöller diagrams. The prescribed time- and 
space-dependent velocity u% (upper left panel). The non-linear advection 
equation integrated numerically with a leap-frog scheme (upper right panel) 
and with the semi-Lagrangian scheme (lower panel). The Courant number 
varies over the interval —1.75 < u < 2.25, which results in a clearly unstable 
solution in the leap-frog case but a nice and smooth solution using the 
semi-Lagrangian framework. 


we can determine the point A’ and compute the interpolated value of 
F in this position: 


Petes -0=p) Fe, (11.18) 


where q = int(2u) and 8 = 2u — q. Here the factor “2” is due to the 
two time steps separating F~} from F/’*?. 
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n-1 


j-q-1 -q j-p-1 j-p 


Figure 11.4. Schematic representing the semi-Lagrangian scheme with 
variable-velocity advection. The “real” (blue curve) and approximate (red 
line) trajectories that arrive at grid point 7 at time level n+ 1. Here a and p 
are the interpolation coefficients used to calculate u2’, 6 and q, the 
quantities used for the interpolation of F”?~?. 


The semi-Lagrangian calculation will thus go through the following 
stages: 


A first guess of u}, which could be u? = u7. 
Compute u, œa and p with this u,. 
Compute a new u* = au”_, ,+(l—a)u 


9=p ' 
Iterate over stages 2 and 3 until no significant improvement is 


AUNE 


obtained. 
5. Compute q = int(2u) and 8 = a —q. 
6. Finally use F?** = Fet = BF} + (1-6) FR. 


The semi-Lagrangian formulation above is only valid for one- 
dimensional problems on a regular grids but can easily be extended 
to curvilinear grids as well as two- and three-dimensional flows. The 
interpolations presented here have all been linear for didactic reasons. 
The disadvantage is, however, that the linear interpolation creates too 
much diffusion of the tracers. The semi-Lagrangian formulations, which 
are employed in circulation models today, for this reason use higher- 
order schemes such as cubic interpolations. The clear advantage of 
the semi-Lagrangian scheme over the Eulerian ones is that one can 
use much larger time steps At. The accuracy of the semi-Lagrangian 
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scheme discusssed in the present chapter is only of first order and will 
decrease with longer time steps. The time step At is often chosen to 
be 5 to 10 times larger for a semi-Lagrangian scheme than for an Eule- 
rian one. A drawback with the semi-Lagrangian schemes has been that 
they were not as mass conserving as the Eulerian ones, which has been 
rectified in the more recent formulations. When applying the semi- 
Lagrangian method to the shallow water equations with forcing, dissi- 
pation and bottom topography some additional difficulties will arise. 


12. Model Coordinates 


In order to present some 3D modelling in the next chapter we will here 
show the different types of vertical coordinates that are used for oceanic 
and atmospheric circulation models. 


12.1 Oceanic vertical coordinates 


The most common vertical coordinate systems used in ocean circula- 
tion models are presented in Figure 12.3. They are of basically three 
types: z-coordinates (constant depth), terrain-following sigma coordi- 
nates and isopycnic coordinates with density layers. The first oceanic 
general circulation model, developed by Bryan and Cox (1967), used 
fixed z-coordinates with a rigid lid as illustrated by the top right panel 
of Figure 12.1. The rigid-lid approximation was replaced (Killworth 
et al., 1991) by treating the fast barotropic mode separately, viz. intro- 
ducing a free surface. Another improvement of the fixed z-coordinates 
was achieved by Pacanowski and Gnanadesikan (1998) by adjusting 
the thickness of the deepest layer in order to match the total depth 
as illustrated by the middle left panel of Figure 12.3. Adcroft and 
Campin (2004) introduced a time dependence in what is known as 
z*-coordinates, where all layers were adjusted to the free-surface ele- 
vation variations like an accordion (middle right panel of Figure 12.3). 
Terrain-following vertical coordinates were suggested by Phillips (1957) 
for atmospheric forecasting models 


12.1.1 Fixed-depth coordinates 


A simple example of a vertical discretisation is the one of the continuity 
equation, which is used in many Ocean General Circulation Models 
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Figure 12.1. Different ocean-model vertical coordinates. 


(OGCMs), based on the B-grid with fixed-depth levels as shown in 
Figure 12.2. The continuity equation 


Ou ðv Ow 


PN le = 12.1 
Ox = Oy = Oz ( ) 
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Figure 12.2. Finite-difference boxes for the B-grid and the C-grid with 
fixed-depth level coordinates. 


can be discretised on a B-grid as 


(Ui jk F Ui,j-1,k) = (Wi-1j,k + Ui-1,j-1,k) 
2Ax 


Wijk = Wign—1 — Az 


(Ui j,k T Ui_-1,5,k) a (Vi j-1,k + Ui—1,j—1,k) 
2Ay , 


(12.2) 


or on a C-grid as 


Uigk — Ui-1,5,6 | Vig, — Vi,j—1,k 
Wijk = Wi,j,k-1 a-( = ee ee . (12.3) 
T y 


Equation (12.2) is integrated from the bottom upwards with the bound- 
ary condition W;,j,x=0 = 0 at the bottom of the ocean. The interpreta- 
tion of this equation is that the sum of all the volume fluxes in or out of 
the grid box is zero. An alternative way to derive the vertical velocity 
is therefore to consider that the sum of the volume transport through 
the six grid-box walls must be zero due to the incompressibility. This 
sum for the C-grid box is hence 


(Ui jk = thagu) AYAZ F (Vije — Vi j-1,x) ATAZ 


(12.4) 
+ (Wijk = Wi jk—1) AzAy = 0, 


which becomes identical to Equation (12.3) by a reformulation. Note 
that we have for simplicity used a k that decreases with depth in order 
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to have the upward positive direction as k increases. In most OGCMs, 
however, k will decrease with depth, with the surface layer counting as 
k=l 


12.1.2 Variable-depth coordinates 

The layer thickness is in most of today’s OGCMs a function of both 
space and time as shown in Figure 12.1. On a C-grid the mass trans- 
ports through the eastern, northern and upper faces, respectively, of 
the i, j, k grid box at time step n are given by 


U? i,k = Piz, kU, KAY j Azi, k? (12.5) 
Vik = Pins pe ATI, Goze gs (12.6) 
Wik = Pig kKWi jk ATi GAY j, (12.7) 


where the unit is kg/s. The continuity equation, which expresses con- 
servation of mass, states that 


dp _ Apu) , Her) , ew) 
ðt Ox Oy Oz 


=0. (12.8) 


Integrating Equation (12.8) over a finite grid box of volume ArAyAz 
we obtain 


ƏM; jr 
at 


where M; jẹ is the mass of the grid box. The rate of mass change of the 


+U, j,k — Ui 5,6 + Vije — Vi,j-1,4 Wize — Wi jr- = 9, (12.9) 


grid box 0M, ;,,/Ot can be either due to compression in a compressible 
GCM or to grid-box volume change, which in a GCM is generally due 
to the time dependence of the vertical resolution so that the thicknesses 
of model layers vary in time. 
The mass of the grid box is 
Min = Pij k ATi j AYi jAk (12.10) 


wJ, 


The vertical mass transport through the top of the grid box is obtained 
by discretising Equation (12.9) between two stored time levels: 


n —_ n n n n 
Wis = Wijk- s- (vz je T Uire t Vije Vie 


ont a i (12.11) 
Pipe aie Pi AZ, 


2At 


i ' 
ahh ANT, Ans) 
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which is computed by integration from the bottom and upwards with 
the bottom boundary condition W;,;. = 0. The vertical velocity w can 
then be deduced from Equation 12.7. 

In many OGCMs, the fluid is considered to be incompressible and 
thus the density is taken to be constant in the equation above. The 
vertical volume transport through the top of the grid box then becomes 


n Aa n n n n n 
Wak = Wijk- 7” k = UF" 1 i,k a Vivi E Vij-1,k 


N , (12.12) 
Az T — Az, 
l 3J; IR 3J; anu) , 


where U, V and W are the volume transports in the unit m?/s. 


12.2 Atmospheric vertical coordinates 


Instead of depth/height as vertical coordinate in our system of 
equations it is possible to use other quantities. The density varies with 
latitude and height/depth which makes the equations sometimes less 
easy to use than an alternative system which employs other quantities, 
such as pressure, sigma or potential temperature for the atmosphere 
and density or sigma for the ocean, as the vertical coordinate. These 
coordinates may facilitate solving the complete equations of motion. 


12.2.1 Generalised vertical coordinates 


We can derive a system of equations for a generalised vertical coor- 
dinate ¢, which is assumed to be related to the height/depth by a 
single-valued monotonic function. When we transform the vertical coor- 
dinate a variable u(x, y, z,t) becomes a(x, y, C(x, y, z,t),t). The hori- 
zontal coordinates remain the same. Let s represent x, y or t. From 
Figure 12.3 we see that 


C-A B-A,C-BAz 


As As | Az As’ piada) 
so that 
Oa Oa Oa Oz 
(a), o) (3). on 
where 
ða Oadz (12.15) 


AG Oz aÇ 
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Figure 12.3. Schematic showing the vertical coordinate transformation. 


or 
ða að 
LS 12.1 
Oz OCOz ( 6) 
Substituting Equation (12.16) in Equation (12.14), we obtain 
Oa Oa ða OC [Oz 
=] = l 12.1 
(a). (Z) +5255 (5). rato) 


From this relationship, we can obtain an equation for a horizontal 
gradient of the scalar a in Å coordinates: 


Oa OC 


= = 12.1 
V:a = V.a + De gz Ve (12.18) 


and for the horizontal divergence of a vector V: 


>- OV a 


iV =V,- reg 2V,z. 12.19 
Ve V Vz e əz Cf ( ) 


The total derivative of a(x, y, Ç, t) becomes 


da 
ac 


12.20 
Di (12.20) 


Da (0a 
Ot 


) BV Ve 
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Figure 12.4. Illustration of pressure, sigma and hybrid vertical coordinates 
in an atmospheric model including terrain. The three types of coordinates 
are superimposed on each other in the lower right panel. 


The three basic (pressure, sigma and hybrid) atmospheric verti- 
cal model coordinates will now be examined and are shown in 
Figure 12.4. 


12.2.2 Pressure coordinates 


Pressure or isobaric coordinates can be used in the atmosphere when 
the hydrostatic approximation (cf. Section 13.1) is applied. In pressure 
coordinates, where 0p/0¢ = 1, the total derivative, Equation (12.20 ), 
is given by 

da Oa 


> Oa 
pare aes ‘ sts 12.21 
TE apy SRE ( ) 
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where the vertical velocity in pressure coordinates is w = dp/dt. The 
continuity equation (12.1) can now be written as 
> Ow 
Vp: V +=— =0. 12.22 
Do (12.22) 
From this expression w can be computed by integrating from the top 
of the atmosphere, where w = 0, downwards over the pressure layers k, 
which on a C-grid yields 


Wijk = Wi,j,e—1 — ADE (= oo Pe a). (12.23) 
This is known as the “the kinematic method” to compute the verti- 
cal motion in the atmosphere. This is basically the same method as is 
used in ocean models expressed by Equation (12.3) and has the clear 
advantage of being mass conserving. The drawback is, however, that 
this motion tends to be very “noisy” since it depends on the diver- 
gence of the wind and hence on its weak ageostrophic component. In 
particular, this becomes a problem when the horizontal velocity field is 
based on observations, as is the case in numerical weather predictions. 
The “adiabatic method” may then be employed instead by using the 
thermodynamic energy equation based on the the first law of thermo- 
dynamics, see e.g. Holton and Hakim (2013). This law, if expressed in 
the isobaric system in the absence of diabatic heating or cooling, is 


+uz—+4v Hw =0. (12.24) 


By defining the the static stability parameter for the isobaric system 
as 


T (12.25) 


where C, is the specific heat at constant pressure, we can deduce an 
expression for the vertical motion: 


1 (OT f OT a OT 
W; 5 — . 
oF (ðt . Ox 7 Oy 


(12.26) 


There are other methods to calculate the vertical atmospheric motion 
such as the “vorticity method” (Sawyer, 1949), which as the names indi- 
cates, is based on a vorticity equation. Most atmospheric GCMs today, 
however, use terrain-following vertical coordinates, described below. 
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Figure 12.5. Two different vertical resolutions of the hybrid-coordinate 
model used at the European Centre for Medium-Range Weather Forecasts 
(ECMWF). 


12.2.3 Atmospheric sigma coordinates 


The sigma-coordinate system defines the origin at the ground level 
of the model. The surfaces in the sigma-coordinate system follow the 
model terrain and are steeply sloped in the regions where terrain itself 
is strongly inclined. The sigma-coordinate system defines the vertical 
position in the atmosphere as a ratio of the pressure difference between 
the location in question and the top of the domain to that of the pres- 
sure difference between the origin and the top of the domain. Because 
it is pressure-based and normalised, it is mathematically easy to cast 
the governing equations of the atmosphere in a relatively simple form. 
The sigma coordinate is hence o = p/ps, where ps(x,y, z,t) is the 
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pressure at the surface of the Earth. The boundary values are conse- 
quently ø = 0 at the top of the atmosphere where p = 0 and o = 1 
at the surface of the Earth. An illustration of the difference between 
pressure and sigma coordinates is presented in the upper right panel of 
Figure 12.4. 


12.2.4 Hybrid coordinates 


The hybrid coordinate system has the properties of sigma coordinates 
in the lower atmosphere and those of pressure coordinates in the strato- 
sphere. Following Simmons and Burridge (1981) the atmosphere is 
divided into Ne layers, which are defined by the pressures at the inter- 
faces between them, these pressures given by 


Pr+i/2 = Ak+1/2 + On41/2 Ps (12.27) 


for k = 0,1,.., Nies, with k = 0 at the top of the atmosphere and k = 
Nies at the surface of the Earth. The a;,41/2 and b;41/2 are constants, 
the values of which effectively define the vertical coordinate and pg is 
the surface pressure. The values a, and b, for a 60-layer model are 
shown in Figure 12.6. 


a (Pa) b (Pa pa) 
a= 0.00000 b= 0.00000 
a = 10.00000 b = 0.00000 
a = 20.00000 b = 0.00000 
a = 28.21708 b = 0.00000 


a = 38.42530 b = 0.00000 


a= 7.36774 b = 0.99402 
a= 3.68387 b = 0.99582 
a= 0.00000 b = 0.99763 


J a= 0.00000 b = 0.99881 


=== MODEL SURFACE i =I a= 0.00000 b= 1.00000 


Figure 12.6. Model levels for the 60 layers shown in Figure 12.5, with 
corresponding a, and by. The middle of each layer in red and the interfaces 
in black. 


Model Coordinates 127 


The dependent variables, viz. the zonal wind u, the meridional wind 
v, the temperature T and the specific humidity q, are defined in the 
middle of the layers, where the pressure is defined as 


1 
Pk = z (Pr-1/2 + Pr41/2) (12.28) 


for k = 1,2,..,Miey. The vertical coordinate is 7 = n(p,pg) and 
has the boundary value 7(0,p;) = 0 at the top of the atmosphere 
and 7(ps,ps) = 1 at the surface of the Earth. Two different vertical 
resolutions with hybrid coordinates are shown in Figure 12.5. 

The vertical mass transport between the layers in a hydrostatic 
AGCM can be deduced from Equation (12.11) so that 


n = n n n n n 
Wijk = We ha — (va. a Oe i,j,k + ijk Vga ih 


m , (12.29) 
Aptt} — Apts} 
ij, LIR Ng. Ay. . 
2gAt Pig Ais J> 


where the hydrostatic approximation has been used. 


12.2.5 Isentropic coordinates 


Isentropic vertical coordinates use the potential temperature, which is 


defined as 
p (Rg/Cp) 
=T (=) , (12.30) 


where R, is the specific gas constant, C, the specific heat at constant 
pressure and po the reference pressure, where 0 = T. 

Isentropic vertical coordinates are convenient when the motion is 
adiabatic since the potential temperature 0 of an air parcel is then con- 
served. In the absence of diabatic processes and mixing, air flows along 
the 0 surfaces, which then act as “material” surfaces. A Lagrangian or 
quasi-Lagrangian vertical coordinate is one that moves with the fluid. 
The isentropic coordinate system therefore qualifies as such. This is its 
main advantage since that “vertical” motion is very weak if the flow 
is quasi-adiabatic, which reduces finite-difference errors in areas such 
as fronts. There are, however, two main disadvantages with isentropic 
coordinates: isentropic surfaces intersect with the ground in contrast 
to sigma coordinates and only statically stable solutions are allowed, 
since the vertical coordinate has to vary monotonically with height. 
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12.3 Structured and unstructured grids 


The finite difference schemes mainly dealt with in this book are applied 
on structured grids such as those illustrated in Figure 12.7. The draw- 
back of the Cartesian grids is that they do not change in space. Oceanic 
and atmospheric circulation models do not have Cartesian grids apart 
from some academic ones with pedagogical aims, cf. the present book. 
A typical grid used for practical purposes will at least have a latitude 
dependence of Az, taking into account that the distances between lon- 
gitudes decrease with increasing latitude. These grids are known as 
curvilinear, but are still orthogonal, i.e. preserve right angles between 
the two coordinates at every point of the grid. Figure 12.8 provides 
an example of this. The region bounded by two adjacent segments 
of one of the curvilinear coordinates and two adjacent segments of 
the other curvilinear coordinate will be transformable to a rectangle. 
An orthogonal curvilinear coordinate system permits the design of a 
grid system with the “north pole(s)” of the coordinate system shifted 
to a terrestrial location. Figure 12.8 shows a tri-polar grid, with the 
two “north poles” shifted to Siberia and the wastelands of northern 
Canada. 

Unstructured grids, in contrast to the structured grids used with 
finite differences, do not require regular connectivity between the grid 
cells. The resolution can hence vary in space, with e.g higher reso- 
lution in coastal regions or narrow straits. Here, we will present a 
brief description of the general principles underlying the finite ele- 
ment and finite volume methods, which both employ unstructured 
grids. 


Figure 12.7. Three examples of structured grids. From left to right: 3D 
Cartesian, 3D rectilinear, 2D curvilinear. Note the “pole problem” of the 
curvilinear grid. 
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Champ en couleur 0); Min= 2.29, Max= 4.75, Int= 0,30. 


700 260 3.20 340 4,40 5.00 5.60 6.20 680 7.40 200 


Figure 12.8. The orthogonal curvilinear ORCA12 ocean grid for the NEMO 
model, which is tripolar with two “north poles” in order to avoid the north 
pole being an ocean point. The colour scale indicates the grid size in km. 


Figure 12.9. Finite elements for an ocean general circulation model 


12.3.1 Finite element method 

The main advantage of the finite element method is that one can locally 
increase the horizontal resolution, which is particularly important for 
coastal regions, cf. Figure 12.9. Disadvantages are, however, that the 
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CF L-criterion will be set by the smallest element, hereby imposing a 
severe limitation on the time step to be used in the integration. Fur- 
thermore data on an unstructured grid can be somewhat problematic 
to analyse. 

The essence of the finite element method (FEM) is recognised by 
considering various ways of representing a function u(x, t) on an interval 
0<a< L. In the finite-difference method the function is defined only 
on a set of grid points; i.e. u;(t) = u(x;,t) is defined for a set of 
zj, denoted nodes, but there is no explicit information about how the 
function behaves between these grid points. 

In the finite element method, the function is defined in terms of 
a finite set of piecewise linear basis functions e,(x) as illustrated in 
Figure 12.10. The variable u(x,t) is assumed to vary linearly between 
the nodes with a piecewise linear fit. The function u(x,t) is then rep- 
resented by the sum 


N 


u(x,t) => uj(é)e;(a), (12.31) 


j=0 


where the grid is defined as x; = 7Ax and Ax = L/N. The basis func- 
tions e; (x) are local, i.e. they are non-zero only on a small-sub-interval. 
Here the u; are the coefficients of the basis functions and u(x) is defined 
everywhere. The 1D linear advection equation employed in Chapter 4 
to illustrate the use of finite differencing, will here be used to describe 
the idea underlying FEM: 


ae M 
ot Ox 


0. (12.32) 


Figure 12.10. Illustration of how the variable u is built up with a 
combination of the piecewise linear basis functions ej. 
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By inserting the piecewise linear representation of Equation (12.31) we 


Ou, ðe; 


where r is the residual. Point collocation, viz. setting r = 0, does 


obtain 


not prove to be useful, but employing the Galerkin technique we may 
impose that 


L 
1 redz = 0 (12.34) 
0 


for i = 0,1,2,...N. Here r can be replaced using Equation (12.33), 
which yields 


ðu; [7 ” ðe; 
2 / exejda + > uj | aede =0. (12.35) 
The following results concerning the basis functions may be derived: 


A 2A 
f (aas; f (Cjpe;)dx =0; f ejdx =; (12.36a) 


de; de; +. 1 de; +. 
—e,dx =0 J A e,dx = +5 | Ge Sdz =0, (12.36b) 


where p is any integer except —1,0,1. These relationships can now be 
used in order to reformulate Equation (12.35) as 


1 duj+1 du; du;—1 Uj+1 — Uj- _ 
; ( itt gaT 4 St |) pe ( SH) 0. (12.37) 


The time derivative of u at the location x; and time t” is, for conve- 
nience, given by 


duj 
G =—2 12.38 
J dt ? ( ) 
which, after substitution into Equation (12.37), yields 
1 n n n u j+1 ` Uj—ı 
7 (Gig AG + G7) = e( 7 Aa ) ; (12.39) 


The right-hand side of this equation is known and consequently it is 
possible to solve this system of simultaneous linear equations for all G7 
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using the methods outlined in Chapter 10. Hereafter Equation (12.38) 
can be integrated in time using e.g. a leap-frog scheme: 
upt! Su + 2AtG?. (12.40) 
The stability of this scheme can be analysed using the von Neumann 
method in the same manner as for the finite-difference schemes. 
A comprehensive overview of oceanic finite element modelling can 


be found in Danilov et al. (2004) and a more fundamental one in Duben 
et al. (2012). 


12.3.2 Finite volume method 


Since the 1980s considerable scientific interest has increasingly been 
directed towards the finite volume method (FVM). Here, instead of 
dealing with variables at specific grid points (as in the case of finite 
differences), we consider their averages over given volumes. The philos- 
ophy behind FVM is that a volume integration of a PDE comprising a 
divergence term converts the latter to a surface integral of fluxes using 
the Gauss theorem. The finite volume method can be applied to both 
structured and unstructured grids. An advantage of FVM over FEM is 
that it conserves the variables better on a coarse grid. 

As in the FEM case, we will use the 1D linear advection equation, 
here reformulated in “flux terms”, to describe the general idea behind 
FVM: 

ðu O(cu) 
Ot Ox 


FVM uses averaged values of the variables over the given cells, rather 


=: (12.41) 


than the values at specific grid points as in the finite-difference method. 
The cell ynt? describes the spatial region between z,_ 3 and z,,1 
and the temporal region between t” and t”*™! (cf. Figure 12.11). The 
spatially averaged value of u over this cell between z,_ 1 and z,, 1 is 


z 1 
m 1 its 


————— ds : 12.42 
pags fue) de (12.42) 


1 
2 
The temporally averaged value of cu over this cell between t”? and t"*! 


1S 
¿n+l 


nth 1 
cu, 2 = af cu(x,,t) dt. (12.43) 
t 
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t= 1 1 i+] £ 
Figure 12.11. AOE A grid for a 1D finite volume model. The 


shaded cell A ta is the spatial region between Ti; and z; Ji and the 


temporal region between t” and t”*!. The “fluxes” in and out of this cell 
are indicated in red. 


Integrating Equation (12.41) over a spatiotemporal cell and dividing 
by AwvAt we obtain 


— zde dt =0. 
Ta [ oe + Asit X k En 
(12.44) 
After integration this becomes 
n+1 
1 Pre apy 1 7 
reas (u u”) de + TG K (ety emg) dt =0. 
e 
(12.45) 
This can be reformulated using Equations (12.42) and (12.43): 
upt up | cunts T CUa o (12.46) 
At Agr l l 
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Once the 1D edge “fluxes” have been calculated using either inter- 
polation or extrapolation of the cell averages, this equation can be 
integrated forward in time: 


upt ut + Ae Cas — a) ; (12.47) 
Although this equation is superficially similar to what may be obtained 
using finite differences, it is important to underline that when physi- 
cal space is multidimensional, the cells can be unstructured and the 
numerical integration clearly differs from that used when applying 
finite-difference methods. An overview of the use of the finite volume 
method in meteorology can be found in Machenhauer et al. (2009). 


13. 3D Modelling 


Ocean General Circulation Models are denoted OGCMs and 
Atmospheric General Circulation Models AGCMs. When coupled to 
each other they are often referred to as atmosphere-ocean coupled gen- 
eral circulation models (AOGCMs). Today’s coupled climate models 
include not only circulation models but also other components such 
as land-surface, aerosols, chemistry, etc. and they are known as Earth 
System Models (ESMs). The core part is still the GCMs of the ocean 
and the atmosphere, which are increasing in complexity as modellers 
improve them. The numerical improvements can be such as higher- 
order numerical schemes, advanced model grids, etc. The improvement 
of the physics can be to replace the hydrostatic equations by the non- 
hydrostatic ones, to change the mixing parameterisations of the unre- 
solved scales, etc. All these improvements are necessary in order to 
continue to undertake more realistic model integrations. 


13.1 Approximations 


A number of approximations and hypotheses are always necessary to 
make in order to simplify the full Navier-Stokes equations. 


The spherical-Earth approximation 


The geopotential surfaces are assumed to be spheres so that gravity, 
coinciding with the local vertical, is a constant. This approximation, 
instead of using the more accurate oblate spheroid, where the Earth 
“flattens” at the poles and “widens” at the Equator as a result of the 
centrifugal force. This approximation is made in all circulation models. 
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The thin-layer approximation 


The thickness of the atmosphere or the ocean is neglected compared to 
the radius of the Earth. This is also sometimes referred to as the shal- 
low atmosphere/ocean approximation, which should not be confused 
with the approximations leading to the “shallow-water equations”. The 
equations without the thin-layer approximation are known as the “deep 
equations” or the “non-hydrostatic deep equations” . 


The vertical Coriolis approximation 


This approximation is sometimes confused with the thin-layer approx- 
imation or the hydrostatic approximation, but is in fact an approxi- 
mation in itself. The 2Qw cos ġ term in the zonal momentum equation 
and the 2Qu cos ¢@ term in the vertical momentum equation are always 
omitted in hydrostatic models, but can be included in non-hydrostatic 
models, which will also be the case for the simple non-hydrostatic model 
presented later in this chapter. 


The incompressibility approximation 


The incompressibility approximation states that Dp/Dt = 0. The mass 
conservation equation expressed by the continuity equation Dp/Dt + 
pv: V = 0 leads to V-V = 0. The 3D divergence of the velocity 
is hence approximated to be zero and the volume is also conserved. 
An accurate expression for density computed by the equation of state 
may, however, be used when computing the pressure if the hydrostatic 
approximation is used. The fluid is therefore often referred to as pseudo- 
incompressible. This approximation is used in many OGCMs but never 
in AGCMs since air is highly compressible. 


The Boussinesq approximation 


Density variations are neglected except when they contribute to the 
buoyancy force, which explains the use of pọ instead of p in the hori- 
zontal momentum equations (13.la) and (13.1b). Since the density is 
still permitted to vary in the vertical momentum equation the approxi- 
mation is therefore more accurately sometimes referred to as the “quasi- 
non-Boussinesq approximation”. “Quasi” since the equation filters out 
the acoustic waves. The Boussinesq approximation is only good if the 
vertical variation of density is small relative to the mean density and 
that the horizontal and temporal variations are small relative to those 
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in the vertical. This is consequently a good approximation for the ocean 
where p ~ 1000 kg/m?. 

It is important to note that this approximation and the incompress- 
ibility approximation both are consequences of the condition that the 
density variations be small compared to the mean density. 


Non-Boussinesq approximations 


The Boussinesq approximation can not be used for a realistic atmo- 
spheric model since the density decreases from 1.2 kg/m? at sea level 
to 0kg/m? at the top of the atmosphere. It is instead possible to use 
other approximations such the anelastic or the pseudo-incompressible 
approximations in order to filter out acoustic waves. 

The anelastic approximations are similar to the Boussinesq approx- 
imation but permit a vertical variation of the mean density so that 
p = po(z)+p' (£, y, z,t), where po(z) is the density satisfying the hydros- 
tic balance and p'(x, y, z,t) is a small perturbation around this balance. 
This eliminates sound waves by assuming that the flow has velocities 
much smaller than the speed of sound and permits a decreasing den- 
sity with height. The anelastic approximations have, however, some 
important limitations since they deform the Rossby modes. 

A more realistic approximation is the pseudo-incompressible approx- 
imation, which accounts for density fluctuations that arise from the 
the equation of state. Density fluctuations associated with perturba- 
tions in the pressure field are neglected. The pseudo-incompressible 
equation is the same as the anelastic continuity equation when the 
mean stratification is adiabatic. When the stability increases, how- 
ever, the pseudo-incompressible approximation gives a more accurate 
result. The pseudo-incompressible approximation is, however, more 
exact when the stratification is stronger and the air more stable. 
The pseudo-incompressible approximation includes the effects of tem- 
perature changes on the density in the mass-conservation equation, 
which is not included in the anelastic approximation. 


The hydrostatic approximation 


The vertical momentum equation is reduced to a balance between 
the vertical pressure gradient and the buoyancy force, which removes 
convective processes from the equations. Convection is instead 
parameterised with an increased vertical diffusion. This approximation 
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is suitable for the large-scale circulation of the ocean and the 
atmosphere. But when the horizontal scale is shorter or on the same 
order as the depth scale this is not accurate anymore. We will, there- 
fore, present a simple hydrostatic model as well as a non-hydrostatic 
model in this chapter. Note that use of pressure coordinates requires 
the hydrostatic approximation, cf. Section 12.2.2. 


The turbulent-closure hypothesis 


The urbulent-closure hypothesis is that the turbulent fluxes (repre- 
senting the effect of small-scale processes on those of larger scales) are 
expressed in terms of large-scale features. In the present book this 
is most often expressed in terms of simple Laplacian diffusion and 
viscosity but can be parameterised in other ways. It is important to 
understand here that a model with higher resolution will not “need” as 
much sub-grid parameterisation as one of lower resolution since the for- 
mer instead will resolve the scales better and “use” the hydrodynamic 
equations. 


13.2 A simple hydrostatic model 


In the present section, we will present a 3D circulation model formu- 
lated as simple as possible, using all the approximations above. A Carte- 
sian grid will therefore be used, which is suitable for the ocean since 
depth is used as a vertical coordinate, this in contrast to atmospheric 
models employing pressure-dependent vertical coordinates. The numer- 
ical core of our 3D model will be close to that of a GCM. It should, 
however, be emphasised here that we do not recommend coding exactly 
these equations since their use would be very limited. The model is here 
discretised only for didactic reasons. Nearly all OGCMs are based on 
some type of curvilinear coordinates and depth-dependent layer thick- 
nesses, which makes the discretised equations less transparent. We have 
therefore discretised the equations for a rectangular Cartesian C-grid 
as illustrated by Figures 13.1, 13.2 and 13.3. 

Our simplified 3D-model is thus on a rectangular ocean domain 
with a flat bottom and furthermore makes use of the hydrostatic 
and Boussinesq approximations: The equations of motion that will be 
used are 
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Figure 13.1. Horizontal view (longitude-latitude) of a possible rectangular 
model grid with land represented as yellow grid boxes. Only the non-zero 
variables are shown. I and J are the total number of grid boxes in the zonal 
and meridional direction, respectively. 


Ou Ou Ou Ou 1 Op Ou 
— E| A F”, 
e a ay ae ne 
(13.1a) 
Ov Ov Ov Ov 1 Op v 
sao l A A ar 
a ae ey on oo uV? v + vaat 
(13.1b) 
i (13. 1c) 
ðu Ov Ow 
aE le eee 13.1 
Oo a oe eas) 
OT T ð ð OT 


o T T 
tug tus twa = KaViT + Kyo, +Q, (13.1e) 
y 


2 
ue gee sige tu% = = KyViz StR s 3 


= p(S,T,p). (13.1g) 


(13.1f) 
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Figure 13.2. Zonal-vertical view of the model grid with land/sea floor 
represnted as yellow grid boxes. Only the non-zero variables are shown. K is 
the total number of vertical depth layers. Note that the vertical index k is 
increasing with depth and is hence in the of opposite direction of that of the 


z-coordinate. 
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Figure 13.3. C-grid with points for the zonal velocity u, meridional velocity 


v, water or air column height h and vorticity €. 
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Here A and K are the viscosity and diffusion coefficients, respec- 
tively. The subscripts H and V pertain to horizontal and vertical 
processes, respectively. F” and F” represent the wind forcing, Q the 
thermal forcing at the sea surface. These seven equations for the ocean 
circulation correspond those established by Bjerknes (1904) for the 
atmosphere. The major differences are that the equation for humid- 
ity has been replaced by one for salinity and that air is compressible, 
whereas water has been taken to be incompressible. To prepare for 
discretisation it is convenient to rewrite the equations as 


Ou ðu OE ru 7 
Ov ðv OF By 
= HA Ay +F”, 13.2 
J Eu aF Dy aVqu + vat (13.2b) 
Tagi (13.2c) 
Oz 
Ow ðu Ov 
= 13.2d 
Oz Ox Oy’ ( ) 
OT O(juT) awT) awT) 5 OT 
= K T+ K 
at a ay. ae eg T 
(13.2e) 
Os O(uS) wS) A(wS) A Os 
= KyViS+K 13.2f 
at i a o ee a 
p = p(S,T,p), (13.28) 
where the absolute vorticity has been defined as 
ðv Ou 
See 1 
Cau (13.3) 
and the energy function as 
=P, (u? +07). (13.4) 
Po 2 


A possible discretisation of these equations with centred finite 
difference is 


urtit = A i 1 + 2At 
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(13.5a) 
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a f- 4 [ i,j,k (ur in + U e) + oF ik (U iik + ti) na] 


n n 
Ui j,k—1 — Vi, j,k+1 


wt: n Tt: nm 
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8Az 
Erik — Eijk 
Ay i 
(13.5b) 
k—1 
Dijk = D IPi jkt AZ + Pij AZ |2 + Pij aij» (13.5c) 
kf=1 
Wi,j,k-1 = Wi,j,k Veg Tiers _ Vag Vasa (13.5d) 


Ax Ay , 


where the absolute vorticity is located between the corners of the 
T-boxes as illustrated by Figures 13.1 and 13.3: 


= Vitigk — Vigk UWi,jtisn — Ui,j,k 
spa fs 3 13.6 
Ei j,k f Ax A y ( ) 


The fluxes U and V are defined at the same points as the velocity 
components u and v: 


1 

Ui j,k = Uijikg (hi jk hizijk) » (13.7) 
1 

Visa = tages (hije + hi jtik). (13.8) 


2 


The grid-cell thickness is constant at all depths except at the surface, 
where the sea-surface elevation h is taken into account: 


hijr = Az +m; fork=1, (13.9) 
hijr = Az for k £1. (13.10) 


The gradient operator will act on the quantity E defined at the same 
locations as h: 
—Pijk , lL |1y o 


1 
E; jk = Pa ] 919 (ix ] U aje) T 5 Cavemen . (13.11) 


13.3 The tracer equation 


The tracer equation describes the rate of change of a tracer such as 
e.g. potential temperature or salt in the ocean, water vapour in the 
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Figure 13.4. Horizontal (left) and vertical (right) views of the tracer 
equation applied on a C-grid. The blue arrows illustrate the tracer flux 

Ui jk = tijk (Ti j,k + Tii, j,k) AyAz and W; j,k = Wij (Tije +Tijk+) 
Ardy. 


atmosphere or any tracer that is advected and diffused in the ocean or 
atmosphere. With a simple parameterisation of the diffusion the tracer 
equation can be expressed as 


oT OT 
a +V.VT = KyV?, T+ Kvo ae (13.12) 


where Ky and Ky are the horizontal and vertical diffusion coefficients 
and Q a possible source term such as the heat flux between the atmo- 
sphere and the ocean. In the case of incompressibility, as postulated 
for our ocean model, the continuity equation is 
> Ou Ov. Ow 
V-V = Fota = 13.13 
Ox Oy Oz ( ) 
The tracer equation (13.12) can now be rewritten by incorporating the 
continuity equation (13.13): 


oT 


P PT 
atv: (V7) eT ke +0 (13.14) 
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13.3.1 Discretisation on a Cartesian grid 

This tracer equation will now be discretised on the C-grid illustrated in 
Figure 13.4. Let us start with the term V- (V7) in Equation (13.14), 
which expresses the divergence of the tracer flux. The discretised ver- 
sion of this term is the sum of all the tracer transports in and out of a 
grid box divided by the volume of the grid box. The tracer flux across 
the grid wall where u; jẹ is located becomes 


1 m m 

OF = Wijk (Trk + Taigin) AyAz, 
1 n m 

Vije = Viik (Tie + Ti ie) AZAZ, (13.15) 
1 n n 

Wir = Wijk (Tir + Tr ac) Ardy, 


which leads to 


(Ui 5% = Ui—1,j,k t V; j,k = V; j—1,k +F W., j k-1 — Wijk) 
AxAyAz l 


V- (Vr) N 

(13.16) 

With a centred leap-frog time scheme and the diffusion terms evaluated 

at time step n — 1 as required to ensure stability, the discretised tracer 
equation becomes 


nm+1 __ n—1 
Tige = Liget 
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(13.17) 


13.3.2 Discretisation on an orthogonal curvilinear grid 


A drawback of the discretised tracer equation above is that it requires 
Cartesian grids, which do not change in space. Oceanic and atmo- 
spheric circulation models do not have grids of this type apart from 
some academic ones used in courses on numerical methods. 
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It is therefore advantageous to apply finite differentiation directly to 
the diffusive tracer fluxes. Let U, V,W now instead be the sum of the 
advective and diffusive fluxes so that 


1 Tr) — Te] 
n hyn n n itljk — tijk 
er ar. (The + Than) Ka Le | AY, g AZ 
=l n—1 
1 TT 
n = n n n ijt 1k i,j,k 
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-1 n—-1 
1 TT — T”; 
n = n n n i,j,k ijik+1 
Wee = eee, (Tije + Theta) — Kve Az AB, jAYi 7. 
k 


(13.18) 


Note that we have written Ky; with an added index k, which indicates 
that we can permit the vertical diffusion coefficient to vary with depth. 
It can be assumed to be either constant, or a function of the local 


ð du\* 
R= 75, (=) (13.19) 


or computed from a turbulent closure model using either the TKE or 
KPP formulation. 
The discretised tracer equation now becomes 


Richardson number 


n+1 __ m—1 
Tit = Tea 
n n n n nm n 
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(13.20) 


Note the additional horizontal indices on the grid lengths Ax; ; and 
Ay;,;, this in order to conform to a curvilinear grid such as the one in 
Figure 12.8. The vertical grid thickness Az, will, however, only have 
a single vertical index since it solely varies in the vertical with the 
exception of the bottom grid boxes that might vary horizontally in 
order to fit an exact depth of the ocean. 


13.4 Non-hydrostatic modelling 


In a similar way as in the previous simple hydrostatic model we will here 
present an example of a simple non-hydrostatic model. This is based 
on the incompressible Boussinesq approximation of the equations of 
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motion in z-coordinates. We still have the same seven variables (u, v, 
w, p, p, S, T), which are computed with the same set of seven equations 
but without the hydrostatic approximation and with the vertical 
Coriolis terms included. 


1 
> sa > F Ge, (13.21a) 
o 10 
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o 10 
S = Sa na ey (13.21¢) 
o o o 
oe 4 y + 5: =0, (13.214) 
OT 
zr = Or (13.21e) 
OS 
p = p(S,T,p), (13.21g) 
where 
G, = u "5 w + 2Q(vsinġ — w cos ġ) + F”, (13.22a) 
ene use og w 2Qusind + F”, (13.22b) 
G,= use v wot + 2Qucosd + T+ F”, (13.22c) 
0 
oT OT OT 2 oT 
Qr = Ur "By Uz T KuVal + Re t+ @a, (13.22d) 
Os Os OS | ag 
Qs = ua, "a wa KuVuS + Kya. (13.22e) 
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© 43.214) ef 2 143,218) $ © 13 016) (13.23) 
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and by using the continuity equation we obtain 
V2p = pV - G, (13.24) 


where G = (G,,G,,G.). This elliptic equation for pressure replaces 
the vertical integration of the hydrostatic equation (13.1b). 
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A difference when solving the governing equations with the 
non-hydrostatic terms included are that the vertical velocity is solved 
directly from the vertical momentum equation instead of integrating 
the continuity equation vertically. Furthermore the pressure is obtained 
by using an equation derived from the continuity equation together with 
the three momentum equations. 

The fundamental difference between an oceanic non-hydrostatic 
model as the one presented here and an atmospheric one is that in 
the atmosphere one does not apply the Boussinesq approximation and 
that the fluid is compressible. The acoustic waves are therefore included 
since they are slow enough to be resolvable explicitly in the horizon- 
tal directions, whereas an implicit treatment in the vertical direction 
is sufficient to make the solution stable. Non-hydrostatic atmospheric 
models are essentially hyperbolic in pressure in contrast to the oceanic 
ones that are elliptic. 


14. Spectral Methods 


In some atmospheric general circulation models (AGCMs), the 
horizontal spatial representation of scalar dynamic and thermodynamic 
fields is based on truncated series of spherical harmonic functions. The 
nature of the underlying two-dimensional horizontal physical grid, also 
known as a transform grid, is tightly coupled to the parameters of the 
spherical harmonic expansion itself. 

The numerical integration methods discussed so far are based on 
a discrete representation of the data on a grid of points encompass- 
ing the space over which a prediction of the variables is desired. Local 
time derivatives of the quantities to be predicted are determined by 
expressing the horizontal and vertical advection terms, sources etc. 
in finite-difference form. Finally, the time extrapolation is achieved 
by one of many possible algorithms, e.g. the leap-frog scheme. The 
finite-difference technique has a number of associated problems, such as 
truncation errors and linear as well as non-linear instabilities. Despite 
these difficulties, the finite-difference method has been the most prac- 
tical method of generating forecasts numerically from the dynamical 
equations. As mentioned above there is another approach known as the 
spectral method which avoids some of the difficulties cited previously, 
in particular non-linear instability; however, this method is less versa- 
tile and the required computations are comparatively time-consuming. 
In a general sense, the mode of representation of data depends on their 
nature and the shape of the region over which the representation is 
desired. An alternative to depiction on a grid of discrete points is a rep- 
resentation in the form of a series of orthogonal functions. This requires 
the determination of the expansion coefficients of these functions, 
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and is said to be a spectral representation, or a series expansion, 
in wave-number space. When such functions are used, the spatial 
derivatives can be evaluated analytically, eliminating the need for 
approximating them using finite differences. 

Many AGCMs employ spectral techniques for their horizontal dis- 
cretisations. This method has become the one most widely used for inte- 
grating the governing equations over hemispheric or global domains. 
The spectral method is, however, also used in some regional NWP 
models such as HARMONIE-AROME (Bengtsson et al., 2017). 

One of the main virtues of the spectral technique is the combination 
with semi-implicit methods. Grid-point models have to use expensive 
methods to solve the elliptic balance equation required at every time 
step in a semi-implicit model. In spectral models this is instead inex- 
pensive and accurate. Furthermore, semi-Lagrangian techniques avoid 
the calculation of nonlinear advection terms which is expensive in pure 
spectral models. Spectral transform models can compute advection 
terms less expensively. Finally the representation of positive definite 
quantities like moisture is an issue in spectral models. 


Associated Legendre polynomials 


The associated Legendre polynomials are the canonical solutions of the 
general Legendre equation 
d d 


dn (1 z’) (2) = -n(n + 1)P”, (14.1) 


where the indices n and m are referred to as the degree and order 
of the associated Legendre polynomial, respectively. The independent 
variable can be reparameterised in terms of angles, letting x = cos, 
where 0 = m /2— ọ is the colatitude. These polynomials to a low degree 
n and order m are presented in Table 14.1 and in Figure 14.1. 


14.1 Spherical harmonics 


As basis functions global atmospheric models use spherical harmonics 
(Tallqvist, 1905), which are the eigenfunctions of the Laplace equation 
on the sphere: 


1 1 ory 1 = ( Ze] = n(n + 1) 


gey” = n 
a? | cos? OA? cosy Op ons Oy 


nm 
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Table 14.1. The associated Legendre polynomials P?” in x = cos6, where 
6 = 7/2 — ọ is the colatitude. 


m=0 m=1 m= 2 m=3 
=0 1 
n=1 cos 6 =n 
n=2  4(3cos?0-— 1) —3 cos 6 sin 8 —3 sin? 0 


=3 4(5cos*@—3cos0) —3(5cos?@—1)sin@ 15cos@sin?@ —15sinë 0 


The spherical harmonics are products of Fourier series in longitude A 
and associated Legendre polynomials in latitude ¢: 


YPO) = PP (ue, (14.3) 


n 


where u = siny, m is the zonal wavenumber and n is the “total” 
wavenumber in spherical coordinates (as suggested by the Laplace 
equation). 

In the usual application of the method, the basic prognostic variables 
are vorticity, divergence, temperature, a humidity variable, and the log- 
arithm of surface pressure. Their horizontal representation is in terms 
of truncated series of spherical harmonic functions, whose variations 
are described by sines and cosines in the east-west direction and by 
associated Legendre polynomials from north to south. The horizontal 
variation of a variable U is thus given by 


UO¢,t)=S> SO UPOY2 09), (14.4) 


n=0 m=—n 


where the spatial resolution is uniform over the sphere. This has a major 
advantage over finite differences based on a latitude-longitude grid, 
where the convergence of the meridians at the poles requires very small 
time steps. Although there are solutions for this “pole problem” for 
finite differences, the natural approach when solving the pole problem 
for global models is by using spherical harmonics. 

It is becoming increasingly common to use what is known as the 
“triangular” truncation of this expansion (Figure 14.2). This is defined 
by M = N = constant , and yields a uniform resolution over the 
sphere. The symbol “TN” is the usual way of defining the resolution of 
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Figure 14.1. The first 15 associated Legendre polynomials P’”. The first 

number on the curves indicates the order of the polynomial and the second 

the degree. The upper panel shows these polynomials with m = 1,2 for 

n = 1,2,3,4. The lower panel shows the polynomials with m = 1, 2,3, 4 for 

n = 1,2,3,4. Each colour corresponds to a separate order m of the 

polynomial Př. 
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Figure 14.2. The triangular truncation in the (m,n) wavenumber space 
delimits a three-cornered region of spherical harmonic modes indicated by 
blue dots. Modes outside of this triangle are indicated with red dots. 


such a truncation; N being the smallest “total” wave number retained in 
the expansion. Given that the Earth’s radius is a, the smallest resolved 
half-wavelength in any particular direction is ma/N (320 km for T63, 
190 km for T106), although the corresponding lateral variation is of 
larger scale. 

Derivatives of a spectrally represented variable U are given 
analytically: 


N n 
a ee ae 


n=0 m=—n 


and 


OU Se ee OP ces 
ap 2 2 Un P A 


n=0 m=—n 


14.2 The spectral transform method 


In a spectral model a variable, €(A, p) is represented by a truncated 
series of spherical harmonic functions. This can be expressed as 


mae D EOD Pr (O) (14.5) 
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where j is the latitudinal index, P™(p,) the associated Legendre 
polynomials and the expansion coefficients €”(y,;) are obtained by a 
Fourier transform of €(A, p). 

The grid-point values are obtained by the inverse transform 


Poa >) Pry), (14.6) 


n=|m| 


followed by an inverse Fourier transform to obtain (A, y). 

In a spectral model, the explicit time steps and evaluation of the 
horizontal gradients are undertaken in spectral space. The tendencies 
of the equations are, however, evaluated in grid-point space. One ben- 
efit when representing the variables in spectral space is that horizontal 
derivatives are represented in a continuous fashion, i.e. no finite dif- 
ferencing is needed to evaluate gradients. The methods used for the 
spherical harmonic transforms are, however, beyond the scope of this 
book. It is therefore suggested to use an already existing library to 
carry out the spectral transforms. 


14.3 The shallow-water equations on a sphere 


The momentum and mass continuity equations governing the motion 
of a rotating, homogenous, incompressible and hydrostatic fluid can be 
written in vector form as 


a = -fk x V- Vð +vV?V, (14.7) 
(i) =f 
° =-6V-V, (14.8) 


where V = (u,v) is the horizontal velocity vector, ® is the geopotential 
height, f is the Coriolis parameter and v is the horizontal diffusion 
coefficient. Furthermore, 


d ð = 
—=—4V-V 14.9 
dot en 
and the V operator is defined in spherical coordinates as 
1 1 
V= Cl (14.10) 


acosyO0A  aðp’ 


where a is the Earth’s radius, A is the longitude coordinate and ¢ is 
the latitude coordinate. 
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The equations above describe the shallow-water equations in the 
u,v, ® system. In the model, we use another form of these equations. 
By introducing the relative vorticity Ç and horizontal divergence ô, 
the equations can be transformed into the the ¢,6,® system. We 
do not derive these equations here since this procedure is relatively 


straightforward. 
By introducing 
€=k-(VxV), (14.11) 
and 
5=V-V, (14.12) 
one can obtain the following set of equations (with u = sin y) 
ze ~ a(1 - L?) x olen iy h et] 


a5 1 8. 18 ; U24 V? 
Z + ———— 14.14 
oaa a ( +s). re 
a 1 ə La 
Feet He (14.15) 


where 7 = E + f is the absolute vorticity, which includes the Earth’s 
rotation and (U,V) = (u, v) cos y. 


15. Theoretical Exercises 


The aim here is to apply the theory given in this book and to hereby 
provide a better understanding of basic numerical methods. 


15.1 Exercises given in the main body of the text 


These exercises are mainly found at the end of each pertinent chapter. 


15.1.1 Finite differences 
Determine the order of accuracy of the centred discretisations of the 
advection scheme 


n+1 a n-1 


J kj 
za V 


“ Hei (15.1) 


Solution: 
The Taylor series of f(x) that is infinitely differentiable at x is a power 
series: 

o 1 oy 
fla) = Feo) +(@—a) E] ie- r) SF 


L=LYQ 


T=“ 
By defining u(x, to) = u7, the following Taylor series is obtained: 


ee = u(zo + Az, ty) = 


2 
u(x, to) + (ae E Ax — go) & +2(ge+ Ag — go)? y 


| zo zo 
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3 
+2 (sgt Ax — x)? £3 


zo 


2 
= u(zo, to) E£ Ax Gu l + 4(Az)? cs 
ae) 


and analogously 
+ 
u, = u(zo, to + At) = (15.2) 


2 
u(Zo,to) + K + At — K) ale + 6+ At — 16)" ca 


(0 


3 
+5 (Wt At — 6) Ss 


+ Ol(f6 + At — 6)*] 


u 2u 
= u(zo, to) + At Beles + (At)? G 


ot2 to 
ELAH e a O[(At)’]. (15.3) 
The temporal derivative is 
n+l n-1 3 
u u; 1 ðu 2 p Pu j 
= 2At — —(At)r — O|(At 
2At =l arl, * 6! ) ae |, 7 ea 
Ou Š 
= ( + O[(At) ) 
(0) 
and the spatial derivative is 
Ugi — Uza 1 Ou Z $ Pu j4 
IAr c (2a a + g (At) as 7 + O|(Az)? | 
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Adding both of these expressions, the discretised advection equation is 
recovered: 


n+l n—-1 n aj 
Y A M ig IA A Ou 


2At 2AT Ot 


Ou 


Le 
F Ox 
0 


+ Ol(Az)?, (At)?]. 


zo 


From the analytical solution it is known that Ou/Ot + cOu/Ox = 0. 
Therefore, the centred scheme in both time and space applied to the 
advection equation has an accuracy of O[(Az)?, (At)?]. 


15.1.2 Stability Analysis 


1. Consider the leap-frog scheme for the advection equation: 


n+l n-1 n n 
uj =u; Uj+ı — Uj- 
= 0. 15.4 
I oh es) 
Use 
u? = Ae Ary, (15.5) 


J 


and show that the amplification factor is 


At ? 
jf aa a iw) 
Az Ax 


2. Show that for u > 1 in the exercise above, one of the solutions 
of the differential equation will “blow up”, at least for some 
wave lengths. 


Solution: 
Introducing a wave solution into the discretised advection equation 
(15.4), the following expression obtained: 


À _ AT! eikAa — eika 
”—Q. 15.7 
( 2At °° 20a ) “3 oo 


Rearranging the terms and looking for the non-trivial solution u? 4 0, 
it is found that, 


+ 2ic® sin(fkAx)A-—1=0 > 
A = —ipsin(kAx) + y1 — (wsin(kAz))’, 


where u = cAt/Az is the Courant number. 
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We will compute the norm of A; there are two cases depending on the 
magnitude of u: 


if |u| <1, JA? = AAt= (—ia + V1 — a?) - (Gat V1 — a?) = 1, 
if |u| > 1, A is purely imaginary. Therefore, 


|A|? = 2? sin? (kAzx) (15.8) 


—17- 2p sin(kAc) / (1 sin(kAx))” — 1. 


In the “worst” case kAgx = Z, |A|? = 2u? — 1 2u/p? — 1 and 


Q? 
since |u| > 1, one of the roots is larger than one, and the 


solution will blow up. 


3. Discretise the advection equation with Euler-forward schemes 
in both time and space. Show that for c > 0 (backward scheme), 
the amplitude of the solutions will grow in time (these thus 
being unstable). But for c < 0 (forward scheme) the amplitude 
will decrease in time, i.e. the solution is stable. 


Solution: 
The discretised advection equation using Euler-forward schemes in both 
time and space: 


2 =Q. (15.9) 


Introducing a wave solution into the discretised equation yields 


A-1 eas] 
no. 15.1 
( 7G +c a )w 0 (15.10) 


Rearranging terms, an expression for À is obtained: 


t): (15.11) 
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In order to analyse the stability of the scheme, the squared norm of A 
is calculated: 


AP = Lae (eter = 1)] - [1 — e$t (ere 1)] 
= 1 — u (e^t —1)—p(e**4* — 1) + u? (1+1 etse — € tk”) 
= 1 — p[2cos(kAzx) — 2] + x? [2 — 2cos(kAz)| 
= 1 — [2 cos(kAx) — 2] (u + p?) 
= 1 + 4sin? (*42) (u + p?) = JAP. 


The stability of the scheme is determined partly by the sign of the 
Courant number p: 


if u > 0, then |A|? > 1 and the solution is unstable, 


if u < 0, then we have two different cases. If u > —1 the 
solution is stable. However, if u < —1 the scheme becomes 
unstable. The stability criterion in this case becomes 


A 
—] < Rass <0. (15.12) 
Az 


4. Undertake a stability analysis of the following discretisation of 
the advection equation: 


uth — tum tu) u u” 
j 2\ jt1 j-1 j+1 j—1 
=0. 15.13 
At a Oke eo 


Solution: 
Introducing the wave solution into the discretised advection equation 
yields 

\— A ae e~ than) ikAw —ikAw 


ë —e 
At 


H =Q: 15.14 
c A ur =0 (15.14) 


Rearranging the terms and looking for the nontrivial solution u? 4 0, 
it is found that, 


A 
A = cos(kAzx) — io sin(kAz). (15.15) 
x 
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15.1.3 Accuracy of the numerical phase speed 


Derive the numerical phase speed 


i xi [cât . 
Cp = gaz 5n E sin(kAv)| : (15.16) 


Solution: 
Start by considering the advection equation with centred schemes in 
both time and space: 


n+1 — pt] 


J J 
i " ~—« he 


= ii, (15.17) 


A wave solution of the form u? = uge**4*~Cp"4 is introduced into 
Equation (15.17). 
The following relationships are to be employed: 


upt! zZ uge VAn Eprat) a tiga VODAT Se ONAT) = e RCD Aty? 
aaa _ perry = tye RCD eh GAr-Opnât) _ eh DAty? 
n å — ik((j+1)Az—-CpnAât) __ ikAg pik(jAx—-CpnAt) _ pikAz, n 
Ujit = Uge = Uge e =e Us) 
u? = uge AP Game = uge Or ene Set = eraeun, 
We obtain the following expression: 
—ikCpåt ikCpAât ikAa —ikAa 
E lai D e —e 
+c u; =0. 15.18 
( 2At 2A7 j ca) 


The nontrivial solution is considered. Applying the trigonometrical 


identity e° — e~** = 2isin(a) yields 


—2tsin(kCpAt) | 2ésin(kAr) _ 
a boa =O. (15.19) 
Rearranging terms, 
At 
sin(kCpAt) = A, sin(kAe), (15.20) 


and applying the arcsin function to both sides of this equality, the final 
solution is obtained: 


At 
kCpAt = arcsin Fe sin(k] ‘ (15.21) 
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15.1.4 Diffusion and friction terms 
1. Undertake a stability analysis for the finite-difference version of 
the Rayleigh-friction equation with the right-hand side of 
Equation (8.3) taken at time step n. 
2. Same as in step 1 but the right-hand side taken at time step 


n-1. 

3. Same as in step 1 but the right-hand side taken at time step 
n+. 
Solution: 


Consider the following equation: 


n+l n-1 
j Uj 


2At 


= yurt”, (15.22) 


where a can be either 0, —1 or 1. Insertion of a wave solution yields 


A)! 


After taking into account the nontrivial solution and rearranging terms, 
it is found that 


A? + 2Atyà t — 1 =0. (15.24) 
Three cases are considered: 


Case 1: a = 0 
AX? +2Atyà—1=0, (15.25) 


with the roots 


A12 = —Aty roa leg (Aty)? + 1. (15.26) 


One of these roots is characterised by |A| > 1. The scheme is hence 
unconditionally unstable. 


Case 2: a = —1 
X +2Aty—1=0, (15.27) 


which has the roots 


à = £1 — 2At. (15.28) 
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The scheme is stable for 0 < Aty < 1. However, for 1/2 < Aty < 1 
the roots become purely imaginary and the solution will be oscillatory, 
which is unphysical. The scheme is hence conditionally stable. 


Case 3: a= 1 
A? + 2Aty\? — 1 =0, (15.29) 
where the roots satisfy 


1 
1+ 2Aty’ 


2 


(15.30) 


Both At and y are positive, and hence 0 < A? < 1, and the scheme 
is thus unconditionally stable. 


4. Calculate the stability criterion for 


Ou Oru 


using the following finite-difference scheme: 


f= ge U a (15.32) 


Estimate an upper limit for At when 
A = 10° m?/s, Ax = 400 km (large-scale horizontal diffusion), 
A=1m?/s, Ax = 10m (vertical diffusion in a boundary layer). 


Solution: 
A wave solution is introduced into Equation (15.32), yielding 


\—1 kA __ 9 4 e-ikAs 
E n | u? =0. (15.33) 


At (Ax)? 


The nontrivial solution is considered. Applying the trigonometrical 
identity e’* + e~** = 2cos(a) yields 


At B At .,(kAr 
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In the “worst” case sin? (kAr/2) = 1. The stable non-oscillatory 
solution is given by the following condition: 
(Ax)? 
0< At < ~——. 15.35 
ae (15.35) 


The results for the two cases are 


A = 10° m?/s, Az = 400km > At™* = 11h. 
A=1m?/s, Ax =10m > At™* = 25s. 


5. The diffusion equation can be integrated using the 
Crank-Nicholson scheme: 


jt+1 | j+1 


(Az? (Az? 


Gea Ale 
At 2 


ae ETa La +a 


Examine the stability of this scheme! 
Solution: 
Introducing wave solutions into equation (15.36) yields 


E —1 A = — 24+ eiar eikAr 2 —) } 


At 2 FA 


(Ar) Ar) IEN: 


J 


The nontrivial solution is considered by applying the trigonometrical 
identity e° + e™* = 2 cos(a), and it is found that 


À i + 2AtA sin? (=) = 1 — 2AtA sin’ (=) . (15:36) 


In simplified form: 


_ 1-2AtAsin? (kAx/2) 
~ 14 2AtAsin? (kAz/2) 


zi: (15.37) 


Even in the “worst” case sin? (kAz/2) = 1, À is always smaller than 
one. The scheme is thus unconditionally stable. However, the scheme 
is implicit. 


15.2 Additional theoretical exercises 


In contrast to the previous exercises, the students will need solve the 
following exercises by themselves. 
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15.2.1 Leap-frog scheme 
Study the leap-frog scheme for the advection equation 


1. Apply the leap-frog scheme (centered scheme with 2nd order 
accuracy in space and time). 

2. Derive the stability criterion. 

3. Discuss the computational mode and how it can be avoided. 


15.2.2 Upwind scheme 


Examine the upwind scheme for the advection equation. 


1. Apply the upwind scheme (uncentered with 1st order accuracy 
in space and time). 
2. Derive the stability condition. 


15.2.3 Euler-forward scheme 


Examine the Euler-forward scheme for the diffusion equation. 


1. Derive the stability criterion for the Euler-forward scheme 
applied to the diffusion equation: 


n+l n n — n n 
Uy Uj _ q% 2u} + uj- 


At (Az)? 


2. Discuss how the time step should be chosen. Why is it good to 
choose a smaller time step than that given by the stability 
criterion? Hint: Study how the amplification factor depends on 
the wavelengths, especially those of highest wavenumbers. 

3. Do the oscillations that may appear, due to this scheme, 
represent a numerical mode? 


15.2.4 Staggered vs. unstaggered grid 


The following code lines are taken from a larger Fortran code that 
interpolates the temperature at a given position (described by a real 
pseudo-index xu). The original code works for an Arakawa C-grid. 
For simplicity we will consider the 1D case where this Arakawa grid 
is reduced to a 1D staggered grid. The xu is defined using u-points 
as a reference while the temperature is defined at an h-point. The 
square brackets enclosing dashed lines should be filled in with pertinent 
instructions! 
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SUBROUTINE interp(xu,T) 

! computes temperature at the position 
! given by xu by interpolation 
IMPLICIT NONE 

INTEGER :: ip, im 


REAL 
REAL 


xu, ax, Tint 
DIMENSION(:), INTENTCIN) ta T 


[ --- ] 


ip = 
im = 
ax = 
Tint 


NINT(xu) + 1 

NINT (xu) 

REAL (ip) - xu 

= T(im)*ax + TCip)*(1-ax) 


[ --- ] 
RETURN 
END SUBROUTINE 


Consider the case xu = 10.4 and the temperature values: T(9) = 5, 
T(10) = 10, T(11) = 20. 


1, 


Compute the interpolated temperature Tint at xu using the 
code above. 

The result obtained for the interpolated temperature Tint is 
wrong, why? “Fix” the code (explain the changes) and 
recalculate Tint. 

Would the original interpolation scheme work correctly (for any 
xu value) if an unstaggered grid was used instead? 


15.2.5 Order of accuracy 


Consider the continuity equation of the 1D shallow-water equations: 


Oh | Ou 


ao (15.38) 


Discretise this equation with centred finite differences on an 
unstaggered grid. 

Derive the order of accuracy of the two finite-difference schemes. 
Repeat the previous tasks using instead a staggered grid. 


15.2.6 Nonrotating 2D shallow-water equations 


Consider the nonrotating 2D shallow-water equations. 


1, 


Discretise these equation on an Arakawa B-grid and apply the 
leap-frog scheme in time. 


168 


2. 


3. 
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Compute the stability criterion for these discretisations. For 
simplicity consider the case Ax = Ay. 
Repeat steps 1 and 2 for the Arakawa C-erid. 


15.2.7 Laplace equation 


Consider the Laplace equation V?® = 0 


I; 
2. 


Discretise this equation in space. 

Set Ax = Ay and write the simplest iterative scheme (Jacobi 
iteration). 

Set up a 4 x 4 grid, with 6 = 1 at the boundaries and iterate 3 
times with a starting state of 6 = 0. 

Repeat the exercise using the faster Gauss-Seidel scheme. 
Repeat the exercise using the even faster SOR scheme. 


15.2.8 Semi-implicit scheme 


Consider the nonrotating and y-independent shallow-water equations: 


1. 


2. 


Ou Oh 

Ov 

=o (15.40) 
Oh Ou 


Discretise these equations on a C-grid using a semi-implicit 
scheme centred at time level n. Terms containing spatial partial 
derivatives must be evaluated as the average at time levels 

(n — 1) and (n +1). 

Undertake a stability analysis of these discretised equations. 


16. GFD Computer Exercises 


The aim of this chapter is to apply our previously gained skills in 
numerical methods to some fundamental problems in geophysical fluid 
dynamics (GFD). 


16.1 Advection and diffusion equations 


The leap-frog, the Euler-forward and the upwind schemes will be stud- 
ied and applied to the advection and diffusion equations. The advection 
equation is 


Ou Ou 
— —=0 16.1 
on on? ot) 
and the diffusion equation is 
Ou 8u 
— — A =0. 16.2 
a op t te) 


The solution interval will in all cases be taken as 0 < x < 1, with the 
periodic boundary condition u(0,t) = u(1,t). 


Relative error 


The relative error E” as a function of time level n is defined as 


F ao am 2N 1/2 
BE" = (Za =u) ) , (16.3) 


D 


where %;’ is the discretised form of the analytical solution and I is the 
number of grid points in the spatial domain. 
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16.1.1 Advection equation 


Examine the the advection equation using a simple numerical model. 
Write a program that can solve the advection equation (set c = 1) with 
two different schemes: the leap-frog scheme and the upwind scheme. 


a) Cosine wave: 


i) Model set-up: Run the program with the resolution 
Az = 0.02 and the Courant numbers 0.9, 1.0 and 1.1. Use a 
cosine wave as initial condition: 


u(x,t = 0) = cos(2rzx). 


ii) Solve the problem using the leap-frog scheme as well as the 
upwind schemes separately. Initialise the leap-frog scheme 
with a single Euler-forward step. 

iii) Plot the results obtained with both schemes and the 
analytic solution in the same figure. 

iv) Analyse the phase error and the amplitude error for both 
schemes. Do the results agree with the theory? (Note: Cp 
can be calculated analytically as the system has a single 
known value for k). 

v) Show how the relative error develops in time (consider 
running the code for as long times as t ~ 104). 

b) Cosine pulse: 


i) Model set-up: Run the program with the resolution 
Az = 0.01 and the Courant number 0.9. Use a cosine pulse 
as initial condition: 


2 


eo { + 4 cos (107 (x — 0.5)) for 0.4 < x < 0.6, 


0 elsewhere. 


ii) Solve the problem using a) the leap-frog scheme and b) the 
upwind scheme. Initialise the leap-frog scheme with a single 
Euler-forward step. 

iii) Try to identify the computational mode. 


16.1.2 Diffusion equation 
Examine the diffusion equation using a simple numerical model. Write a 


program that can solve the diffusion equation using the Euler-forward 
scheme. Set Ax = 0.05, A, = 1, and try with three different time 
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steps: on determined by the critical value of the von Neumann number 
v = AAt/ (Ax)? permitted by the stability criterion, one slightly larger 
and one half of the maximum possible value. 


Run the program with two different initial conditions: 


i) Rectangular pulse: 


a 


i 1 for0< z< 0.5, 
0 for05<a< 1. 


ii) Spike: u = 1 at the single grid point x = 0.5 and u = 0 at all 
other grid points. 


16.2 1D shallow-water model 


Here we will examine the difference between implementing a staggered 

and an unstaggered grid applied to the 1D shallow-water equations as 

well as how to incorporate open boundary conditions in the model. 
The shallow-water equations in the one-dimensional case are 


Ou Oh 
Z ey a 16.4 
ðh Ou 


These equations can be combined into a single expression for u and h: 


2 2 2 2 
o'u gr oh g?r 


oe a — op 97 ae ~ 9: (16,8) 


The general d’Alembert solution to the wave equation is given by two 
functions travelling in opposite directions: 


ulz, t) = u(x — ct) + u(x + ct), (16.7) 
h(a, t) = hi (a — ct) + ho(x4+ ct), (16.8) 


where c = \/gH is the propagation speed of the waves. The solution 
interval to be analysed here is 0 < x < 1, with the periodic boundary 
condition u(0, t) = u(1,t). 

Examine one-dimensional gravity waves, which can be described 
by the shallow-water equations. Write a program that solves the 1D 
shallow-water equations given above. 
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a) Unstaggered grid: 

i) Model set-up: Discretise using the leap-frog scheme on an 
unstaggered grid. Set H = g = 1 (thus c= /gH = 1) to 
simplify the system. Run the program with Az = 0.025 and 
the Courant number 0.9. 

ii) Initialise the leap-frog scheme with a single Euler-forward 
step, and use the following initial conditions: 


5+ 4 cos (107r(x — 0.5)) for 0.4 < x < 0.6, 


0 elsewhere. 


h(x,t = 0) = i 
u(x,t =0)=0. 


iii) Do the results agree with the analytical solution of the 

problem? 
b) Staggered grid: 

i) Model set-up: Rewrite the program using a staggered grid, 
with h and u defined at different points. 

ii) Run the program and try to find the stability limit setting 
the time step. Then choose a time step 10 % smaller than 
that in the stability limit, and employ the same initial 
condition as in case a). 

iii) Run three different simulations: one with the resolution 
given in case a), a second one with half that resolution 
(finer grid) and a third one with double the case-a) 
resolution (coarser grid). 

Select one of the variables (either h or u) and compare the 
accuracy of the solution with the previous solution obtained 
on an A-grid. Discuss the difference between the three 
results and the result from case a). Note that Az is the 
distance between two h-points (and also two u-points). 

c) Open boundary conditions: 

In order to gain an understanding of how to deal with open 

boundary conditions requiring a “sponge” (relaxation) zone we 


iv 


ar 


will start by examining the 1D case. Consider a relaxation zone 
defined in the region 0 < x < L,. The relaxation function y has 
the following values at the boundaries of the sponge zone: 


y(L,) = 0, (16.9) 
¥(0) = 1. (16.10) 
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Here we shall consider two types of relaxation functions: 
Linear function: 


yz) =1-—. (16.11) 


Trigonometric function: 


(2) = ; Í + cos (=) (16.12) 


i) Model set-up: Consider the 1D shallow-water model on a 
domain —L < x < Lona staggered grid. Set c= L = 1 and 
Az = 0.025. Consider as initial condition a centrally located 
cosine-shaped “bump”: 


1,1 107a L Ë 
4+ 1 cos ( ) In — 4, SS 


(h(x, 0), u(a,0)) = (ho, uo) f i T 


elsewhere 


Here ho = uo = 1. Impose a solid boundary at x = L and an 
open boundary at x = —L by adding a relaxation zone 
adjacent to the western boundary. The relaxation can be 
parameterised in the following way: 


ualz) = i — O u(x) for —L,<a<-—L (16.13) 


where T is a relaxation time-scale (s~') and uga is the 
damping value. Construct the function so that y = 0 in the 
interior of the system, and y — 1 as you approach the 
boundary. The function may be either trigonometric or 
linear. Here it is sufficient to set 7 = At for simplicity. 

ii) Analyse the effects of the length of the relaxation zone L,., 
choose values within the interval 0 < L, < L. Select a 
proper magnitude to study the effect of the open boundary. 

iii) Discuss the results and the differences you observe when 
you impose the trigonometric and the linear relaxation 
function. 


16.3 2D shallow-water model 


Here the aim is to solve the linearised 2D shallow-water equations 
using some standard numerical methods and to study atmospheric 
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and oceanic processes. These equations, including rotation and physical 


parameterisations, are 


Ou Oh 

a = fu- 95, +F”, (16.14) 
0 Oh 

S sujang + FY, (16.15) 
Oh ðu Ov 

a -H & a x) (16.16) 


Here the physical parameterisations F® and F” may for instance be of 


the horizontal viscosity and/or the drag forcing. 


Below you will find the general structure of a Fortran-coded shallow- 
water model. The dashed lines should be filled in with relevant Fortran 


instructions! 


PROGRAM structure_of_code 


General structure of a shallow-water model in 


Fortran with staggered grid, 


rotation, diffusion 


It is good and common practice to always write a 


| 
| 
! & relaxation schemes. 
| 
| 


few lines about 


Think about: 
Explain each 


Do not write 


! 
| 
! Use common notation. 
] 
| 


what the code does and how. 


parameter in words and units. 


too much in one line. 


Physical constants: 


REAL*4, PARAMETER f = ? ! Coriolis parameter [s-1] 

REAL*4, PARAMETER g =? ! Gravity [m s-2] 

! Model parameters: 

REAL*4, PARAMETER H =? ! Mean depth [m] 

REAL*4, PARAMETER mu = ? ! Diffusion coeff. [m2 s 
=i] 

! Grid 

INTEGER*4, PARAMETER NX = ? ! Number of i- 
points 

INTEGER*4, PARAMETER NY = ? ! Number of j- 
points 

INTEGER*4, PARAMETER NT = ? ! Number of time 
steps 

REAL*4, PARAMETER dx = ? ! Zonal grid spacing [m] 
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REAL*4, PARAMETER :: dy = ? ! Merid. grid spacing 
] 

REAL*4, PARAMETER :: dt = ? ! Time step [s] 
! Save data to file 
CHARACTER*200 :: outFile = ? ! Name of output file 
! Work variables 
INTEGER*4 :: ic, ip, im, jc, jp, jm, nc, nm, np 
REAL*4 :: du, dv, dh 
! Data matrices 
REAL*4, DIMENSION(NX,NY,NT) :: u, v, h 


! Main loop 
! You might want to indent the code inside the loop 


to 
! make it easier to see where the loop starts and 
ends. 
DO ne HsH=S>= 
np nc+1 
nm = nc-í1 
DO jo=------- 
jp = jc+1 
jm = jc-1 
DO i¢==—=-=>> 
ip = ic+1 
im = ic-1l 
! Reset 
du = 0. 
dv = 0. 
dh = 0. 
t Coriolis 
du = du + ~------------------- 
dv = dv + eheee 
! Sea surface height gradient 
du = du + Senenin 
dv = dv + ------------------- 
! Continuity 
dh = dh + ~-------------7------ 
! Diffusion 
du = du + ------3------------- 
dv = dv + Aneen 
! Relaxation 
du = du + ------------------- 
dv = dv + ------------------- 


! New time step 
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WIC jeit) = s--=-sssesH=ese-SesHessee5 
vic, jce, nit) = -s--css nen 
ulic jerit) = -------- 5-5 n nnn nnn 
! Asselin filtering 
WiC, jc, NC) = asesan ae 
Vie, je,nC) < —oa-sSs-SSo eS sos ess Se Sas 
hic, je; ne) = sese-sSsseeeqsS-oS esse Sass 
! Update time index 

ENDDO 

! Don’t forget to take care of boundary 


conditions, 


! zonally periodic, no-slip, or other. 


! Save the data to file 


END PROGRAM 


The following tasks should be undertaken: 


i) 


iii) 


Model set-up: Here we develop a 2D shallow-water model that 
includes rotation on an Arakawa C-grid. Set F” and F” to 
zero. For simplicity prescribe your domain as —1 < x,y < 1 
where Az = Ay = 0.025, g = 1, and H = 1. Choose a proper 
Courant number and introduce an Asselin filter. Integrate the 
first time step using an Euler-forward scheme and the 
subsequent time steps with a leap-frog scheme and keep in 
mind that periodic boundaries will affect the computation of 
the Coriolis terms. 

Consider as initial conditions: 


u(x, y,0) =0, (16.17) 
v(x, y,0) = 0, (16.18) 
h(x, y, 0) = hoe iw -0w (16.19) 


where L = 1, Lw = L/7 and ho = 1. Impose solid boundaries 
everywhere. 

Run the model and analyse the results for the non-rotating 
case (f = 0). 

Calculate the potential and kinetic energies of the system. 
Study the time evolution of both their magnitudes and check 


GFD Computer Exercises 177 


to what extent the total energy is preserved. Does the Asselin 
filtering affect this result? 

v) Repeat the same tasks for the case of constant rotation, 
f =10-+s~!. Interpret the results physically. 


Note that for each of the subsequent tasks given in this chapter you 
will need to programme new parts in the model code to deal with 
different physical processes such as waves and circulation in the ocean 
and atmosphere. 


16.4 Geostrophic adjustment 


The purpose here is to understand what controls the evolution of a 
disturbance initially at rest on a rotating plane, a process commonly 
known as geostrophic adjustment. To study this phenomenon, you will 
use the 2D shallow-water model developed in the previous section. 


Ou Oh 

OE fov Iar’ (16.20) 
Ov ðh 

— — 16.21 
Oh ðu Ov 
—+Do|=—+=]=0. 16.22 
att (e+e! i me) 


Realistic values of g, Do and f should be used in order to gain under- 
standing of real atmospheric and oceanographic processes. The follow- 
ing should be programmed in the code: 


i) Use parameters that are appropriate to the mid-latitude North 
Atlantic when running the experiments, if nothing else is given, 
viz. L = 5-108 m, Dp = 4000 m, g = 9.81 ms~?, fo = 1074 s7?. 

ii) Prescribe the number of grid cells Nx and Ny in the x- and 
y-directions respectively. 

iii) Let Az, Ay be determined by L and Nx, Ny . 
iv) Let At be determined by a Courant number, phase speed and 
max(Az, Ay). 

v) Let Nr be determined by Tmax, which is the simulation period. 

vi) Construct the main time loop so that the model does not store 
the fields every time step. 

vii) Analyse the effect of the boundary condition by prescribing 
different widths of the sponge zone. 
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16.4.1 Geostrophic adjustment for a step-function disturbance 


Here you will study the geostrophic adjustment for an initial distur- 
bance consisting of a discontinuity in A that is symmetric in y. 


Model set-up: 


i) Implement a disturbance in h, which is described by a step 
function in the zonal direction. A good practice is to 
programme in Fortran the initial condition as a CASE or 
LOGICAL to be able to turn it on or off for other exercises. 

ii) The code already describes a Gaussian disturbance of the 
sea-surface height as an initial condition, do not remove this, 
instead make this as one of the options for the initial condition. 

iii) Use periodic boundaries by attaching the northern boundary 
to the southern. 

iv) Use open boundaries in east and west. Make sure that the 
sponge zone is wide enough. 


Exercises: 


i) Derive the final steady state of the sea-surface height and the 

velocities starting from the linearised shallow-water equations. 

ii) Run the model until the system only varies very little in time 
(steady state). 

iii) Verify your model results with the results you obtained 
analytically. 

iv) Tweak the parameters. Which parameters are most important? 

v) Discuss the results obtained by tweaking the parameters. 


16.4.2 Geostrophic adjustment for a Gaussian disturbance 


Here you will study the geostrophic adjustment of an initial disturbance 
described by a Gaussian perturbation of the sea-surface height. You are 
supposed to study the various energetics of the system. 


Model set-up: 


i) Use a Gaussian disturbance instead of the step function 
employed in the previous subsection. The disturbance should 
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be smaller than the domain, but make sure it is larger than the 
resolution. 

ii) Programme the potential, kinetic and total energies of the 
system. Use the Fortran command LOGICAL to be able to turn 
the computation on or off. Remember that you have to save 
these fields to be able to study them. 

iii) Use open boundaries (sponge) everywhere. 


Exercise 1: 


i) Write down the analytical expression for the total kinetic and 
potential energies. 
ii) Run the model with two different depths; one shallow 
Do ~ 500 m and the other deep Do ~ 10000 m. 
iii) Study how the energies of the system (kinetic, potential and 
total) vary in time for the two cases with different depths. 
iv) Give a physical explanation of the results. 


Exercise 2: 


The system can be considered to be in steady state when the energy 
of the system changes very little with time. 


i) Run the model long enough so that the system reaches to a 
steady state for the depths used in Exercise 1. 
ii) Run the model again using different depths and sizes of the 
disturbance. 
iii) How do the results change? Try to explain the changes and 
compare with theory. 


16.5 Kelvin wave 


The aim of this assignment is to gain an understanding of the influence 
solid boundaries exert on geophysical flows. Although applicable to 
some situations in the atmosphere, the influence of lateral boundaries 
is more important in the world oceans. In this assignment you will use 
a simple shallow-water model with solid boundaries, and look at how 
an initial disturbance evolves, and how a system initially at rest evolves 
as you add forcing. 
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16.5.1 Gaussian disturbance 


Model set-up: 


i) Programme an infinite coast by setting an open boundary on 
the north side of the domain and a solid boundary on the 
south side. 

ii) Impose periodic boundary conditions in the x-direction. 

iii) Programme a new initial condition that describes a 
reasonably-sized Gaussian disturbance of the h-field centered 
in the middle of the coast, so that the partitioned disturbance 
has its maximum value at the wall. 

iv) Introduce a reduced gravity in the system. 


Experiment: 


i) Run the model for two different cases; the North Atlantic 
(ZL =1-10’ m , H = 1000 m) and the Baltic Sea (L = 1- 10° 
m,H=30™m). 

ii) Do any particular kind of waves show up? 


) 
iii) How can you identify these waves? (phase speed, shape, ...). 
iv) Connect your results to the theory. 

v) Run the model again for both cases (North Atlantic and the 
Baltic Sea) using different settings, e.g. changing the 
resolution. 

vi) Look at a cross-section of the wave, and estimate its phase 
speed. 


vii) Compare your results with the theory. 


16.5.2 Equatorial 6-plane 


In this part of the assignment you are supposed to study waves that 
appear on an equatorial 6-plane. Model set-up: 


i) Change the initial condition to a Gaussian disturbance centred 
in the middle of the domain. 

ii) Programme a new Coriolis parameter so that it corresponds to 
an equatorial 8-plane (f = By). Try to take advantage of the 
Fortran command LOGICAL here so that you easily can change 
between a 3-plane or an f-plane. Place the equator (y = 0) 
centrally in the domain. Remember to introduce reduced 
gravity in the code! 
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i) Run the model using values corresponding to the Pacific Ocean 
(L = 2-10" m , H = 1000 m) 
ii) Why do you need a reduced gravity for this simulation? 
iii) Identify the waves in the model. 
iv) Do the waves correspond to theory? 
) 


v) Can these waves be observed in the real ocean? 


16.6 Oceanic Rossby waves 


The aim of this assignment is to study and understand the kind of waves 
that are generated when D(y) or f(y) are “sloping”, viz. D = Do+ay or 
f = fo + By. For this assignment you should use a model with periodic 
and open (sponge) boundary conditions. The initial disturbance should 
be in geostrophic balance to avoid gravity waves. 


Model set-up: 


i) Use an initial disturbance in geostrophic balance. (Hint: you 
have already programmed the initial condition for h in a 
previous exercise. But you need to programme u and v so that 
they are in geostrophic balance.) 

ii) Use periodic boundaries in x by attaching the eastern 
boundary to the western. 

iii) Use open boundaries (sponge) on the southern and northern 
boundaries. 

iv) Programme a sloping a-plane D = Do + ay using LOGICAL (or 
CASE) so that you can choose to have it on or off in the 
simulation. 

v) Programme a new Coriolis parameter f = fo + By. Use the 
Fortran command LOGICAL (or CASE) so that you can choose 
to turn it on or off. 


16.6.1 Rossby waves on a {-plane 


i) Consider a rectangular basin with L = 7- 10° m and 
Do = 4000 m in the mid latitudes (e.g. the North Pacific ). 
Run the model for at least 30 days. 

ii) Start by deriving the phase speed and group velocity for 
Rossby waves in this linear system. (The derivation should not 
be included in the report, only the final solution). 
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iii) Run the model with a G-plane and a constant topography Do. 
iv) Describe and explain the evolution of the system. 
v) Connect the results to theory. 
vi) What kind of waves develop? 
vii) Do they have any distinguishing properties? 


16.6.2 Phase and group velocities 
i) Rerun the model but with different wave numbers. (Hint: To 

change the wavenumber, change the width of the disturbance). 

ii) Compare the obtained phase speed and group velocity to the 
theoretical values. 

iii) Repeat this exercise but keep f constant and vary the 
topography using an a-plane. Discuss the differences between 
the results. 


16.6.3 6 — a compensation 
i) Run the model using a varying topography and f-field 
(D = D, + ay and f = fo + By). 
ii) Calculate the value of a that cancels the 8 effect. What 
happens to the initial disturbance under these circumstances? 


16.7 Atmospheric Rossby waves 


The aim of this assignment is to study and understand the kind of 
waves that are generated when f(y) is “sloping”, viz. f = fo + By. The 
effect of including a zonal mean flow will also be discussed. For this 
assignment you should use a model with periodic and open (sponge) 
boundary conditions. The initial disturbance should be in geostrophic 
balance to avoid gravity waves. 


Model set-up: 


i) Use an initial disturbance in geostrophic balance. (Hint: you 
have already programmed the initial condition for h in a 
previous exercise. But you need to programme u and v so that 
they are in geostrophic balance.) 

ii) Use periodic boundaries in x by attaching the eastern 
boundary to the western. 
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iii) Use open boundaries (sponge) on the southern and northern 
boundaries. 

iv) Introduce a zonal mean flow Uo in the system using LOGICAL 
(or CASE) so that you can choose to have it on or off in the 
simulation. 

v) Programme a new Coriolis parameter describing f = fo + By. 
Use LOGICAL (or CASE) so that you can choose to turn it on or 
off. 


16.7.1 6—plane 


i) Consider a rectangular basin with L = 7- 10° m and H = 4000 
m in the mid latitudes. Run the model for at least 30 days. 

ii) Start by deriving the phase speed and group velocity for 
Rossby waves in this linear system. (The derivation should not 
be included in the report, only the final solution). 

iii) Run the model with a 6-plane and without a mean flow 
(Up = 0 m/s). 

iv) Describe and explain the evolution of the system. 

v) Connect the results to theory. 

vi) 


vill 


What kind of waves develop? 
Do they have any distinguishing properties? 


16.7.2 Phase and group velocities 


i) Run the model as above, but with different wave numbers. 
(Hint: To change the wavenumber, change the width of the 
disturbance). 

ii) Compare the obtained phase speed and group velocity to the 
theoretical values. 


16.7.3 The effect of the zonal mean 


i) Increase the zonal length of the domain to L, = 28-10° m 
while using the same Az as in the previous experiments. 

ii) Run the model using a 3 plane and four different zonal mean 
flows (choose zonal flows in the range 0 < Up < 15 m/s). 
Calculate the group velocity for each case. 

iii) Using the results from ii) and try to find the value of the zonal 
mean flow at which the Rossby wave becomes stationary. Is 
the initial condition preserved? If not, explain why. 
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16.8 Gyre Circulations 


The aim of this assignment is to study the Sverdrup, Stommel and 
Munk theories for barotropic large-scale ocean circulation. To study 
this phenomenon, you will use the 2D shallow-water model including 


friction: 
a = fu gon | AyV?ut+Tu, (16.23) 
- =—fu 15 H ApV?0+Tv, (16.24) 
Oh ðu Ov 
an a(> | =) (16.25) 


where Ay is the horizontal viscosity coefficient and I the Rayleigh- 
friction coefficient. 


Model set-up: 


i) Let the system initially be at rest. 
ii) Use solid boundaries. 
iii) Programme a wind forcing that is constant in x and sinusoidal 
in y. 
iv) Use a -plane 


Do not forget to compute the barotropic stream function in order to 
compare your results with the analytical solution. 


16.8.1 Sverdrup solution 
i) Run the model without any friction and use parameters that 
are appropriate for the North Atlantic. 
) What effect does 6 = Of /Oy have on the system? 
iii) What effect does the wind forcing have? 
iv) Do the results resemble theory? Why/Why not? 
v) Do the results improve if the size of the domain is increased? 


11 


(e.g. rerun the experiment for the North Pacific) 
vi) What are the limitations of the Sverdrup theory? 


16.8.2 Stommel solution 
i) Now run the model using the same ocean basin (either the 
North Atlantic or the North Pacific) with the Rayleigh linear 
friction applied SOLELY to v. 
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ii) Give a physical explanation of this kind of friction. 

iii) Describe and explain the modelled circulation using the 
Stommel theory. 

iv) Rerun the model using a different value of the 
Rayleigh-friction coefficient I. 

v) Describe how the results have changed. 


16.8.3 Munk solution 


i) Now run the model using the same ocean basin (either the 

North Atlantic or the North Pacific) with Laplacian friction. 

ii) Give a physical explanation of this kind of friction. 

iii) Describe and explain the modelled circulation using the Munk 
theory. 

iv) Rerun the model using different values of the horizontal 
viscosity coefficient Aj. 

v) Describe how the results have changed. 

vi) Compare the Stommel and Munk solutions! 
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